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Recent trends of ah initio studies and progress in methodologies for electronic structure 
calculations of strongly correlated electron systems are discussed. The interest for developing 
efficient methods is motivated by recent discoveries and characterizations of strongly correlated 
electron materials and by requirements for understanding mechanisms of intriguing phenomena 
beyond a single-particle picture. A three- stage scheme is developed as renormalized multi-scale 
solvers (RMS) utilizing the hierarchical electronic structure in the energy space. It provides 
us with an ah initio downfolding of the global band structure into low-energy effective models 
followed by low-energy solvers for the models. The RMS method is illustrated with examples of 
several materials. In particular, we overview cases such as dynamics of semiconductors, tran- 
sition metals and its compounds including iron-based superconductors and perovskite oxides, 
and organic conductors of a^-ET type. 
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1. Introduction 

Since the foundation of quantum mechanics, under- 
standing and predicting properties of condensed matter 
from microscopic basis have continuously been a great 
challenge of modern science and technology. Behaviors 
of many electrons primarily determine the diversity and 
rich variety of materials in our environment with poten- 
tial applications for future technology. At the same time, 
many electron systems have been a source of challenges 
of our intelligence on nature, because of their interacting 
and quantum mechanical nature. 

Among all, density functional theory (DFT)^'^ offers 
a standard method for calculating electronic structure of 
real materials. Local density approximation (LDA)^ and 
generalized gradient approximation (GGA) offer practi- 
cal ways, and reasonably predict physical properties in 
a wide range of materials if we consider cases with weak 
electron correlations such as semiconductors. 

However, in materials with strong correlation effects 
such as transition metal compounds, organic conduc- 
tors and rare earth compounds, these methods not only 
lead to a quantitative inaccuracy but also to a qualita- 
tively wrong answer. The most famous example is the 
mother materials of copper oxide superconductors such 
as La2Cu04, where DFT predicts a good metal with a 
half- filled band while La2Cu04 is a typical and good 
Mott insulator with a gap amplitude of about 2 eV.^'^ 

In LDA, the many-body Schrodinger equation is re- 
placed by the Kohn-Sham equation 



1 



^k,j = ^k,j^k,j, (1) 



which contains the external potential ]4xt coming from 
nuclei and the Hartree term of the electron-electron 
Coulomb interaction Vh- The eigenfunction and the 
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eigenvalue with the momentum k and other quantum 
number j such as orbital indices are denoted by il)k^j 
and efe,j, respectively. By solving this single-particle 
Schrodinger equation, the ground state energy and the 
charge density are obtained. Here, the electron corre- 
lation effect is accounted by the exchange correlation 
potential Vy^c- From the Hohenberg-Kohn theorem,^ in 
principle, the solution of the Kohn-Sham equation gives 
the exact ground-state energy of the many-body system 
by a functional of the electron charge density.^ However, 
since we do not know how to treat Vy^c exactly, we resort 
to approximations, for instance by LDA. 

In LDA, the exchange correlation potential is replaced 
by results of a uniform electron gas such as quantum 
Monte Carlo calculations etc.^ Therefore, it becomes a 
good approximation when the electron density does not 
have large spatial variations, which is justified in a good 
metal with electron wavefunctions extended uniformly 
in space. However, in strongly correlated electron sys- 
tems, where electrons become nearly localized and elec- 
tron density fiuctuations are large with large spatial de- 
pendence, the approximation becomes poor. 

Typical strongly correlated electron systems are found 
in transition metal compounds, where the Fermi level Ep 
crosses d bands. Rare earth compounds with /-electron 
bands crossing Ep and organic conductors with p bands 

Ef are also well known correlated electron systems. 
A characteristic and common feature of strongly corre- 
lated electron systems is that their bands crossing the 
Fermi level have narrow bandwidths. The origin of the 
narrow bands is that the spreads of p and / orbitals are 
relatively small as compared to lattice constants, which 
makes the overlap of two orbitals each on the neighboring 
atoms small. The relatively small spreads also make the 
local electron interactions large. Experimentally, these 
compounds are often insulators and "bad metals".^ 

In addition, competitions of tendencies for various or- 
ders and fiuctuations such as magnetic, charge and super- 
conducting orders invalidate mean-field treatments in- 
cluding LDA. Electron correlations have to be treated 
at much higher level of accuracies. This is a grand chal- 
lenge of first-principles calculations for electronic struc- 
ture. Strongly correlated electron systems have attracted 
interest as platforms of possible innovative devices and 
realizing functions and efficiencies beyond the semicon- 
ductor apphcations in the 20th century. Ah initio meth- 
ods hold a key of clarifying basic properties from the 
scientific points of view. 

Strong correlation effects appear not only in typical 
correlated electron materials, but also show up even in 
weakly correlated systems such as semiconductors, if ex- 
citations and dynamics are involved. This typically 
emerges in excitonic effects, where an electron and a 
hole interact strongly with attractive interactions. Dy- 
namical ffuctuations also generate effective interaction 
in a small energy scale such as van-der-Waals interaction 
and dispersive forces, where the force is mediated by dy- 
namical electronic polarizations due to electron correla- 
tion effects. Such dynamical ffuctuations are beyond the 
tractability of LDA. However, these weak forces play es- 
sential roles in solutions, complex systems and biological 



systems. For example, they are crucial in determining 
structures of proteins and DNA. In this article, these dy- 
namical effects on dispersive forces are not discussed in 
detail, though it remains a challenge. 

Methods of electronic structure calculations can be 
classified into two categories. The first one is the den- 
sity functional theory described above, where the ground 
state is obtained only from the charge density. The other 
is the wavefunction method that explicitly seeks for so- 
lutions of many-body wavefunctions. One of the sim- 
plest wavefunction methods is, though not sufficient, the 
well known Hartree-Fock theory. Variational wavefunc- 
tion method and Monte Carlo method based on the 
path integral may also be regarded as wavefunction ap- 
proaches. 

The density functional theory reduces the problem to a 
single-particle one through solving the Kohn-Sham equa- 
tion. Here, the self-consistent equation is reduced to ob- 
taining the charge density and the computational load 
is much smaller. On the other hand, the wavefunction 
method allows more ffexibility of treating the electron 
correlation effects, giving us more information on the 
ground state while it is in general more time consum- 
ing even for the Hartree-Fock level. Within the limited 
computer power, the density functional method has thus 
been used more widely. However, from the incentive for 
treating the electron correlation effects more accurately, 
the serious limitation of the density functional theory 
revealed recently has urged reexaminations of standard 
methods. In fact, the standard density functional theory 
so far offers no ways of systematic improvements on this 
electron correlation problem. 

Electron correlation effects can be taken into account 
by considering the standard many-body perturbation 
theory, if the correlation effects are not too large. When 
they interact each other, one can view that each elec- 
tron is dressed by other electrons from the single-particle 
picture and moves in the cloud of other electrons. The 
dressed electron is called a quasiparticle, where the 
particle-like identity is retained in a gedanken experi- 
ment of switching on the interaction gradually in an adi- 
abatic fashion. The quasiparticle may, however, have a 
mass and a dispersion different from a bare electron. This 
is manifested by the self-energy of electrons S, where the 
pole of the quasiparticle (dispersion) is renormalized as 
uo = e*(A:, u) with e*(/c, u) = e(/c) + E(A:, cj), depending on 
momentum k and frequency uj. In the lowest order per- 
turbation, S is given by GW^ where G is the the electron 
Green's function G(k, cj) obtained from the Kohn-Sham 
Hamiltonian and W is the screened Coulomb interac- 
tion. The screened interaction is calculated based on the 
random phase approximation (RPA). This is called GW 
approximation as we describe details in §2.3.^^^^ This 
theory has succeeded in improving the gap of semicon- 
ductors and insulators as we see in §5.1. 

To circumvent the failure of reproducing band gaps in 
correlated insulators, the LDA+U method has been de- 
veloped. ^^^^ It combines LDA with an artificially intro- 
duced "/7" term which raise the energy of electrons only 
for the unoccupied part to take into account the onsite 
Coulomb interaction in the same spirit as the Hartree 
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Fock approximation. 

Recently, more thorough efforts have been made to 
overcome the difficulty of DFT by combining and uti- 
lizing the flexibility of the wavefunction method for a 
more accurate description of electron correlation effects. 
This hybrid approach allows reducing the heavy com- 
putational task of the wavefunction methods and simul- 
taneously allows accurate solutions by improving wave- 
function methods. In this review, we figure out recent 
studies along this line and discuss achievements as well 
as future perspectives. 

The density functional theory is formulated to give 
ground state energies as a functional of the electron den- 
sity only. By extending this, several attempts have been 
made to represent not by the electron density functional 
but by functionals of more information. For example, one 
attempt is to represent by a functional of the whole elec- 
tron Green's function G(k, uj) depending on the momen- 
tum k and the frequency u. The density functional theory 
can be formulated as a theory to minimize a functional 
r(n, Vxc) of the electron density n(r) and the exchange 
correlation potential Vxc • By extending it and replacing 
by an extremum problem for a functional of the Green's 
function G{k^uj) and its conjugate field X, one can have 
a formalism equivalent to that by the Luttinger-Ward or 
Baym-Kadanoff functional. As we discuss later, the dy- 
namical mean field theory can be regarded as one of these 
attempts. The GW method can also be formulated by the 
Luttinger-Ward formalism. In addition, an attempt for a 
formalism including the two-body correlation functions 
in the functional has also been made.^^ 

The origin of the failure of LDA in treating strongly 
correlated electrons is ascribed to reconstructions of elec- 
tronic states near the Fermi level taking place beyond 
the expectation by LDA. The electrons whose energies 
are far away from the Fermi level are either fully oc- 
cupied or empty and do not have a polarizability even 
in the strongly correlated systems. Since the electronic 
polarizability is large near the Fermi level, electron cor- 
relation effects appear in the energy window around the 
Fermi level in the order of the effective electron inter- 
action (, which is typically several eV). Therefore, when 
the widths of bands near the Fermi level become compa- 
rable or even smaller than the effective interaction, the 
whole band structure of this band may be seriously re- 
constructed. Since this happens in the whole Brillouin 
zone, the reconstruction may happen locally, namely in 
a spatially inhomogeneous fashion. This invalidates the 
applicability of LDA around the Fermi level. In other 
words, the LDA may give an adequate band structure 
in the global energy scale, while it is seriously recon- 
structed near the Fermi level in the range of the effective 
interaction (~ several eV), which constitutes a hierar- 
chy structure in energy. Meanwhile physical properties 
of materials around or below the room temperature are 
determined in this low-energy part of the hierarchy. 

By considering this hierarchy structure, one can de- 
velop a ffist-principles method that starts from the den- 
sity functional theory, and then eliminates the degrees 
of freedom far away from by following the spirit of 
the renormalization group. This method of eliminating 



degrees of freedom and restricting the Hilbert space is 
called the downfolding method. ^^^^ The downfolding 
leaves an effective model represented only by the degrees 
of freedom near E^. Since the effective model contains 
only a small number of bands near the Fermi level, for 
which we call "target band" , it is constructed on the lat- 
tice in real space with this number of retained orbitals in 
the unit cell, as in the multi-band Hubbard-type models 
in Lagrangian forms in general or in Hamiltonian forms if 
the retardation effects caused by the downfolded (elim- 
inated) bands are small. Recently, effective low-energy 
models obtained after the downfolding have extensively 
been employed and solved by accurate low-energy solvers 
to discuss strong correlation effects. 

The whole procedures constitute the three-stage 
scheme of the renormalized multi-scale solvers (RMS). 
We here summarize the present RMS method for the 
electronic structure calculation as the hybrid-type three- 
stage scheme as we illustrate in Fig.l: 

(1) Calculate the global band structure including bands 
far from the Fermi level by relying on DFT such as 
LDA, GGA or GW. 

(2) Perform the renormalization procedure to downfold 
the higher energy degrees of freedom. This yields 
effective models for low-energy degrees of freedom 
near the Fermi level. 

(3) Solve the low-energy effective model by an accurate 
low-energy solver. 

Although we do not describe, a possible iterative proce- 
dure to feed back the solution of the low-energy solver 
into the global structure in the ffist step taken until 
the self-consistent solution is a future issue when the 
low-energy solution seriously modifies the original global 
structure. 

This article is organized as follows: In §2, several points 
useful for the ffist stage of the three-stage RMS scheme in 
calculating global electronic structures are summarized. 
After an elementary remark on DFT in §2.1, we intro- 
duce basically three basis functions developed for solving 
Kohn-Sham equation. The ffist is the plane-wave basis 
suited for sp electron systems, and the second and the 
third are the augmented wave and the muffin-tin orbital, 
respectively, suited for d and / electron systems. In §2.3, 
we review GW method as a method to take into account 
electron correlation effects for the global band structure 
within a perturbative approach. In §2.4, the LDA+U 
method proposed to implement Hartree-Fock level cor- 
rections to LDA is reviewed. In §3, the second stage of 
the three-stage scheme of RMS is introduced for the pur- 
pose of downfolding and eliminating the degrees of free- 
dom far from the Fermi level. Effective low-energy mod- 
els are derived from this downfolding procedure. Tools 
for solving the derived effective models are listed in §4, 
where we review, dynamical mean- field theory in §4.1, 
many-variable variational Monte Carlo method in §4.2, 
and path-integral renormalization group method in §4.3. 
Section 5 describes some applications to various materi- 
als as dynamics of semiconductors, transition metal com- 
pounds, and organic conductors. Section 6 is devoted to 
summary and future scope. 
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Fig. 1. (Color online) Schematic procedure of three-stage scheme for hierarchical electronic structure reviewed in this article. 



2. Global Electronic Structure 

2.1 Density Functional Theory (DFT) 

Understanding properties of matter from first princi- 
ples is a central problem in condensed matter physics. 
The properties are, in principle, described by the many- 
body Hamiltonian, 



H 



ij 



■R/l 



--y 



ZiZje^ 



(2) 



where the first term represents the kinetic energy of elec- 
trons. The second, third and fourth terms are interac- 
tions between electrons, electron-nucleus, and nuclei, re- 
spectively. Electrons are labeled by real space coordinate 
r with suffices with the lower case as i and j, while nuclei 
are denoted by coordinate R with the upper-case suffices 
as / and J. The electronic bare mass and charge are m 
and e, while the atomic number is denoted by Z. The 
spin degrees of freedom, relativistic effects and quantum 
effects of nuclei are neglected for simplicity. The Hamilto- 
nian is solved exactly only in very limited cases, hence de- 
veloping a practical procedure for treating many-electron 
systems has long been an important issue. 

DFT gives an approximate but reasonably accurate 
and practical method for this problem. DFT is based on 
the Hohenberg-Kohn (HK) theorem,^ that asserts: 

Theorem 1) For any many electron systems under the 
influence of an external potential ]4xt(r), the po- 
tential is, apart from a trivial additive constant, a 
unique functional of the one electron density of the 
ground state. 

Theorem 2) For any external potential, there exists a 
total energy functional of one electron density n(r). 



^totN = F[n] + / Fext(r)n(r)dr , 



(3) 



where F[n] is a universal functional of n(r). The 
ground state energy of the many electron system is 
the minimum of £^totN, and associated n(r) is the 
electron density of the ground state. 

The HK theorem was originally proved for systems hav- 
ing non-degenerate ground state. Later on it was ex- 
tended to degenerate cases by Levy.^^ The theorem is 
an exact theory of interacting many electron systems. 
Since the Hamiltonian is determined by the ground state 
electron density, all properties of matter are implicitly 
determined by the density. This gives a justification to 
take the electron density as a basic variable of the theory. 

In DFT, the ground state total energy and density 
are obtained by minimizing the total energy functional 
with respect to n(r). The formulation may be regarded 
as a rigorous extension of the Thomas-Fermi (TF) the- 
24, 25 which the total energy functional is given as 



iTFr 



= T 



TF 



Fext(r)n(r)d^ 



1 / n(r)n(r0^3^^3^/ 



1 f n(r)n(r 

2 J |r-r' 



= l(37r2)2/3 I „(r)5/3d3r . 



(4) 



(5) 



The kinetic energy in the TF theory is approximated as 
the integral of the local part over space, where the local 
part is the mean kinetic energy per electron multiplied 
by the electron density at the position. The TF theory 
was proposed in the 1920's and applied to real materi- 
als. However, the method turned out to be unsatisfactory 
not only quantitatively but also qualitatively: The the- 
ory cannot describe chemical bonds between atoms. It 
was clarified that the error comes mainly from the ap- 
proximation for the kinetic energy. 

Much better results are obtained by replacing the ki- 
netic term with that of the noninteracting electron sys- 
tems. This is nothing but the Hartree theory which was 



J. Phys. Soc. Jpn. 



Invited Review Paper 



Imada and Miyake 5 



developed soon after the TF theory. Inspired by this 
observation, in 1965 Kohn and Sham^ proposed a prac- 
tical procedure for DFT.^^ They introduced an auxiliary 
noninteracting electron system that obeys the following 
single-particle equation (Fig. 2) 

|-^V2 + «eff(r)| 4>j{r) = ej^,{r) . (6) 

The electron density and the kinetic energy of the system 
are 

occ. 

n(r) = ^|V,(r)|^ (7) 

j 

OCC ^ 

T. = j2(^j\-^v'm. (8) 

j 

Coming back to the original interacting system, the func- 
tional F[n] in eq.(3) can be divided as 

F[n] = TM + \j ^y^dVdV + E^M ■ (9) 

The first term is the kinetic energy, but it is for the non- 
interacting system defined in eq.(8), not the true kinetic 
energy of the interacting electrons. The second term is 
the electrostatic energy (Hartree energy). The last term, 
so-called exchange-correlation energy, contains all the re- 
maining contributions including the difference between 
the noninteracting and interacting kinetic energies. Now 
we assume that the ground state electron density of the 
interacting system can be represented as eq.(7). Then, 
the stationary condition for the total energy functional 
eqs.(3) and (9) is satisfied when the self-consistent solu- 
tion of eq.(6), with the effective potential 

is achieved. The set of equations (6), (7) and (10) is called 
Kohn-Sham equation. 

The remaining question is how to determine the 
exchange-correlation energy functional. First of all, the 
exact functional is not known, and trials to improve the 
functional is a hot topic even today. Formally the func- 
tional can be written as 

^xcN = j exc(r; [n])n(v)d?r , (11) 

where exc(r; [n]) is the exchange-correlation energy per 
electron at the position r. In principle, full information 
of the density n, not only the value at r is necessary to 
determine exc- A simple and most widely used approxi- 
mation is the LDA proposed in the Kohn-Sham work.^ 
The LDA approximates exc to be that of a uniform elec- 
tron gas of the density at the position. Namely, the en- 
ergy functional is expressed as follows. 

^xc^^W = / exc(n(r))n(r)dV . (12) 

The explicit formula for exc has been proposed by several 
authors based on a perturbation theory, the RPA,^^ 
or more accurately by the fit^^'^^ to the Ceperley- Alder 
quantum Monte Carlo simulation.^ 



The LDA is by construction exact in the limit of a uni- 
form electron density, whereas the approximation gets 
worse as the spatial variation of the electron density be- 
comes strong. Typically the lattice constant of solids and 
bond lengths between atoms are computed to be within 
the 2-3 % error to experiments. The accuracy of the ion- 
ization energy in molecules and cohesive energy in solids 
is with 10-20 % errors. The high accuracy is partially ra- 
tionalized by the fact that the LDA satisfies a sum rule 
for the exchange-correlation hole.^^ 

An obvious modification of LDA is inclusion of density 
gradient effects Vn(r). However, it turned out that the 
simple low-order gradient expansion does not improve 
but worsen the results in some cases. This suggests the 
spatial variation is so strong in real materials that sim- 
ple gradient expansion does not work. Instead it may 
be a better idea to include the effect of gradient cor- 
rections while keeping various asymptotic behaviors and 
sum rules, such as the one for the exchange-correlation 
hole mentioned above. This is called the generalized gra- 
dient approximation (GGA), 

Eg^M = j e^«^(n,Vn)n(r). (13) 

There are many explicit formula of GGA proposed by 
today. ^^^^ The GGA tends to give more accurate results 
than the LDA in the atomization energy, cohesive energy, 
description of magnetism and so on. 

While the total energy and the electron density are ob- 
tained from the total energy functional, the Kohn-Sham 
equation merely represents a fictitious system which is in- 
troduced to carry out the minimization. One may want 
to regard the eigenvalues {cj} as the orbital energies. 
However, this interpretation is not justified rigorously. 
Physical meaning of the Kohn-Sham energy is known 
only for the highest occupied state. It is proved that —Cj 
for the state is the ionization energy. For other states, a 
similar relation Cj = dEtot/dfj holds, but it is not the 
electron addition/removal energy. Here fj is the occupa- 
tion number of the state j, namely, n(r) = fj [^/^^(r)^. 
Mathematically the Kohn-Sham eigenvalue is a Lagrange 
multiplier corresponding to the orthonormal condition 

In practice, the Kohn-Sham energy is useful infor- 
mation for understanding the electronic properties. The 
overall feature of the electronic structure is captured in 
LDA/GGA, and the Fermi surface is reasonably accurate 
in many materials. However, the band gap of semicon- 
ductors and insulators are underestimated significantly. 
This is the case for almost all materials including weakly 
correlated systems. The low-energy electronic structure 
in strongly correlated materials are often very different 
from measurement, and sometimes qualitatively wrong. 

2.2 Basis Functions for DFT 

The first step of describing the low-energy properties 
of correlated materials is the global electronic structure 
by means of DFT in LDA, GGA or whatever. Various 
numerical techniques have been developed to solve the 
Kohn-Sham equation accurately and efficiently. Below we 
will see a few key ingredients. 
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Fig. 2. (Color online) Schematic representation of density functional theory 



2.2.1 Plane wave 

One of tlie most widely used basis functions is the 
plane wave basis set, where the wavefunction of an elec- 
tron with the momentum k and quantum number n are 
expanded as 



^kn(r) 



G 



G gi(k+G).r 



(14) 



A great advantage of the plane-wave basis is that nu- 
merical accuracy is improved systematically by increas- 
ing the number of plane waves (G points in eq.(14)). 
However, the calculation becomes tremendously heavy if 
a naive plane-wave expansion is adopted, because huge 
number of plane waves is required to express localized 
core electrons. To make calculations feasible, the inter- 
actions between core and valence electrons are replaced 
with a pseudopotential. The pseudopotential eliminates 
explicit treatment of the core electrons from the Kohn- 
Sham equation. The valence orbitals are also modified 
to be smoother than the true (all-electron) ones near 
the core region, which reduces computational cost dras- 
tically. The concept of the pseudopotential dates back 
to the 1930's.^^ It has been developed continuously, 
and non-empirical pseudopotential appeared in the late 
1970's. 

The pseudopotential is constructed in such a way 
that the scattering properties of valence electrons repro- 
duce those of the all-electron calculation accurately. The 
pseudo wavefunction agrees with the all-electron wave- 
function outside a certain radius Tc, whereas for r < Tc, 
the pseudo wavefunction is nodeless and much smoother 
than the all-electron one. The pseudopotential can be 
written in the following form. 



Vps(r) = V\oca\{r) 



Im 



m 



(15) 



The first term is independent of angular momentum / and 
is called local part. The second term is dependent on 
and called non-local part. Since (i) the pseudo wavefunc- 
tions are equal to the all-electron ones at r > Tc and (ii) 
the all-electron potential is independent of /, it follows 
that 5Vi = at r > Tc- Various ways of pseudopotential 
construction have been proposed so far to increase ac- 
curacy, transferability and computational efficiency. ^^^-'^ 
For more technical details, see e.g. Refs.42,43. 



2.2.2 Augmented plane waves and muffin-tin orbital 

While the plane-wave basis is suitable for sp electron 
systems, for systems containing d and / electrons, the 
computational cost becomes heavy because of localized 
nature of the electrons. All electron methods are pow- 
erful in such cases. The augmented plane waves (APW) 
and muffin-tin orbital (MTO) are commonly used basis 
functions. 

In the APW method, space is divided into two parts: 
the muffin-tin region and the interstitial region (Fig. 3). 
Inside the muffin-tin region, the basis function is con- 
structed by solving the Schrodinger equation for the 
spherically symmetrized potential at a particular energy 
e. The solution (j) is connected at the muffin-tin surface 
to the plane wave. Thus, the APW is expressed as 



Xk+G(r;e) 



gi(k+G).r > R) 



^Zma(r;e) ,(r <R) 
(16) 

where a is the index for a muffin-tin. 

The Kohn-Sham eigenvalue is obtained as a solution of 
the secular equation — e^*! = 0, where S is the overlap 
matrix between the APW basis functions. Because the 
APW is implicitly energy dependent, the secular equa- 
tion is a nonlinear equation. It is not a general eigenvalue 
problem, but one has to search for the selfconsistency of 
e numerically, which is computationally demanding. In 
1975, Andersen proposed a linear method to solve this 
problem. The energy dependent (j)ima{^] e) is expanded 
at around a fixed energy e/^, and approximated as 



(r; eia) + B (r; eia) 



(17) 



where (p is the energy derivative of (f). The coefficients 
A and B are determined from a matching condition at 
r = up to the first-derivative. The linearized function 
is called linear augmented plane wave (LAPW). It is the 
most accurate method among electronic structure meth- 
ods available today. 

Muffin-tin orbital (MTO) is another basis function for 
the all-electron method. It is defined by 

{(/>zma(r; e) + cot Sia{K)ji{Kr)Yim{0, (f), 
(r < R) 
ni{Kr)Yim{0,ip), {r>R). 

(18) 

Here ji is the spherical Bessel function, and ni is the 
spherical Hankel (Neumann) function for tz'^ > (/^^ < 
0). 6ia is determined by matching the logarithmic deriva- 
tive at the muffin-tin boundary. The MTO is not the 
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Interstitial region 



Fig. 3. In both the APW and MTO methods, the space is divided 
into the muffin-tin region and the interstitial region. The basis 
function is constructed by approximating the potential to be 
spherically symmetric in the muffin-tin region, and constant in 
the interstitial region. 

eigenfunction of a single muffin-tin potential because of 
the second term for r < R in eq.(18), but it is useful for 
solving the many muffin-tin problem. The second term 
expresses approximately the tails of the MTO's centered 
at other sites. It reduces the number of required basis 
functions compared to other methods. 

The Hnear muffin-tin orbital (LMTO)"^^ is the linear 
method for the MTO, in which the following approxima- 
tions are adopted. Firstly the energy e is fixed. It makes 
the basis function energy independent. Secondly, jiYim 
in eq.(18) is replaced with 

cot Sia ' 

This form is chosen to satisfy the following condition at 
the fixed energy. 

^Xw(r) = . (20) 

Thirdly, the niYim term is replaced with a linear combi- 
nation of (pima in other muffin-tins. 

Further efficiency is achieved by approximating the 
whole space as a sum of muffin- tins. In this atomic sphere 
approximation (ASA), spherical anisotropy of the poten- 
tial inside the muffin-tins is neglected. The muffin tin 
radius is chosen so that the interstitial region becomes 
small, while the overlap between muffin-tins is small as 
well, since both the interstitial region and the overlap 
region is neglected in the ASA. 

Although the ASA is not accurate in materials with 
large empty space or strong anisotropy, the LMTO- 
ASA is a computationally cheap and powerful method 
for closed-pack and localized-electron systems. The ASA 
makes mapping onto lattice models easy. Therefore, the 
LMTO-ASA has played an important role in the develop- 
ment of electronic structure techniques for strongly cor- 
related materials. For example, both the LDA+U and 
the LDA+DMFT methods were developed on top of the 
LMTO-ASA in the beginning and extended to other ba- 
sis sets later. 

To go beyond the linear approximation and improve 
the accuracy, NMTO was developed recently. In the 
LMTO, (j) is computed at a fixed energy, whereas the 
NMTO basis is a linear combination of N such functions 



evaluated at TV different energies. 

2.3 GW approximation 

Many-body perturbation expansion is a traditional 
theory for interacting electron systems. Expansion in the 
Coulomb interaction, v{r) = l/|r|, gives the Hartree- 
Fock (HF) approximation^^ in the lowest order but in 
a self-consistent manner. The HF self-energy is a sum 
of two terms: the static Coulomb interaction (Hartree 
term) and the exchange interaction (Fock term). Assum- 
ing that the one-electron wavefunction is not modified 
by electron addition or removal, it can be shown that 
the eigenvalues of the Hartree- Fock equation are electron 
addition or removal energies (Koopman's theorem). In 
other words, the eigenvalue is equal to the total energy 
difference between the N and ± 1 electron systems. 

The HF theory is a good approximation in finite sys- 
tems, but the accuracy goes down in extended systems. 
In fact, the result is often even worse than the Hartree 
theory in solids. The Hartree theory yields too small 
band gap of insulators. In the HF theory, the exchange 
term pushes down occupied energy levels and widen the 
band gap. The correction is, however, too large, con- 
sequently the HF overestimates the band gap. Another 
well-known drawback in the HF theory is the anomaly 
at the Fermi level. The Fermi velocity in metals di- 
verges, and the density of states vanishes at the Fermi 
level. (Comparing to the DFT-LDA, the HF approxima- 
tion is computationally more demanding because of non- 
local Fock term, and the results are worse in extended 
systems.) These facts suggest that proper treatment of 
screening effects is crucial in solids. A sensible way would 
be the expansion in the series of the screened Coulomb 
interaction. This is the basic idea of the GW approxima- 
tion.^^^^ 

Before the GW method is established, there are a few 
works reported in the 1950's. Quinn and Ferrel studied 
the electron gas, and attempted to include correlation 
effects in the form of the GW approximation, with sev- 
eral other approximations.^^ DuBois did a related work 
on the electron gas in high density region. Baym and 
Kadanoff also mentioned a GW form of the self-energy in 
their paper on the conserving approximation.^^ In 1965, 
Hedin derived an exact closed set of equations for the self- 
energy in which the self-energy was expanded in powers 
of the screened Coulomb interaction. In particular, the 
first term in the expansion yields the GW approximation, 
which can be viewed as the time-dependent Hartree ap- 
proximation for the self-energy. 

In his 1965 paper, Hedin presented a full self-energy 
calculation for the electron gas. Shortly after that, 
Lundqvist did extensive calculations of the electron gas 
self-energy and spectral functions. The calculation is 
so heavy that it took 20 years before the applications to 
real materials were reported. Hybertsen and Louie per- 
formed the GW calculation of semiconductors and found 
that the GW approximation cures the band gap problem 
of DFT.^^ Their seminal work used the plasmon-pole ap- 
proximation, which is a simplification for the frequency 
dependence of the dielectric function. Soon after this, 
Godby, Schlueter and Sham carried out the GW cal- 
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culation without the plasmon-pole approximation, and 
obtained similar results. These calculations were per- 
formed using the pseudopotential methods based on the 
plane- wave basis. However, pseudopotential GW calcu- 
lations for materials containing localized electrons are 
computationally demanding. This difficulty motivated 
all-electron GW calculations. The all-electron GW calcu- 
lations were done by Aryasetiawan in 1990's.^^'^^ With 
the rapid increase in computer performance, GW cal- 
culations can now be performed for systems containing 
more than 50 atoms. 

The basic quantity of the GW approximation is the 
one-particle Green's function defined by 



G{l,2) = i{N\T[i;{l)i^H2)]\N) 



(21) 



where \N) is the ground state of the N electron system, T 
is the time-ordered product, ip and i/j^ are field operators, 
and 1 = (ri,ti) is the composite variable. Starting from 
the equation of motion of the Green's function, Hedin 
derived a set of equations between the Green's function 
G, self-energy E, screened Coulomb interaction W, po- 
larization function P, and vertex function F: 

E(l, 2) = ^ y G(l, 3+)V^(14)r(3, 2, 4)d(34) , (22) 

W(l,2) = ^(l,2) + z J ^(1, 3)P(3, 4)1^(4, 2)d(34), 

(23) 

P(l, ^) = -iJ G(l, 3)r(3, 4, 2)G(4, l+)d(34) , (24) 

r(l,2,3) = (5(1-2)^(1-3) 
(5E(1,2) 



/ 



(5G(4,5) 



G(4, 6)G(7, 5)r(6, 7, 3)d(4567) , 



(25) 



G(l,2) = Go(l,2)+ / Go(l,3)S(3,4)G(4,2)d(34). 



(26) 

A key to solving the equations is the vertex func- 
tion. To solve the integral equation (25), we must know 
STi/SG^ which requires an explicit expression for the self- 
energy in terms of the Green function. But the self-energy 
depends in turn on the vertex as can be seen in (22). In 
the GW approximation, the vertex function is approxi- 
mated as 



r(l,2,3) = (5(l-2)(5(2-3) 



which leads to the following form of the self-energy (, to 
which the name of the approximation owes). 



E^^(l,2) = zG(l,2)W^(l,2) . 



(28) 



The diagrams for eqs.(23) and (28) are illustrated in 
Fig. 4. We note that the Hartree contribution in the self- 
energy shown in the first line of Fig. 9 is not considered, 
because it is already considered in the Green's function of 
Kohn-Sham equation in the LDA level. The polarization 



function is expressed as 

P(l,2) = -iG(l,2)G(2,l+) 



(29) 



The GW approximation is the lowest order expansion of 
the self-energy in the screened Coulomb interaction. Con- 
sidering that the Fock term is i;x(l, 2) = iG(l, 2)v{l — 2), 
GW may be regarded as the screened Hartree-Fock ap- 
proximation using the screened Coulomb interaction in 
the RPA. 



w 



« = = = • 

1 2 



1 2 1 3\t^4 2 



2 = 



Fig. 4. Feynman diagrams included in the GW approximation. 
Solid lines with arrows represent Green's function represented 
by eigenstates of Kohn-Sham equation. Dashed lines represent 
the bare Goulomb interaction v, while double dashed lines are 
for the screened interaction W. 



Most applications to real materials are carried out non 
self-consistently on top of the LDA solution. The starting 
Green's function equivalent to eq.(21) is 



Go(r,r» = ^ 



Mr)rAr') , "-^'^c(r)^*(r') 

uj — Cv — i5 



E 



UJ — Cc -\- iS ' 

(30) 

where tjjy and tjjc are the eigenstates of Kohn-Sham equa- 
tion. The polarizability eq.(29) is then expressed as 



Po(r,r^c^) 



occ unocc /' I 

ee{- 



(r)V^*(r)V^,(rOV^:(rO 



^ T f . (31) 

UJ ^ [Cc- Cy) - id 

from which W is computed. Putting the obtained W 
into eq.(28), E is computed. Representing eq.(28) in the 
frequency domain, the self-energy can be divided into 
two terms. One is the contribution from the poles of the 
Green's function 



(27) EsEx(r,r» = -^tA«(r)tA:(r')W^(r,r';e„-w). (32) 



This is the same form as the Fock potential in the HF 
approximation, but the Coulomb interaction is replaced 
with the screened one. It is called screened exchange 
term. The other term comes from the poles of called 
Coulomb hole term. 



all ^, 

S]coH(r,r» = VV^(r)V'*(r')P / 

Jo 



UJ — Ci — UJ' 



(33) 

where B is the spectral function of W. The meaning of 
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this term is understood by taking the hmit, 

uj-Ci^O . (34) 

In this static Coulomb hole + static exchange (COH- 
SEX) approximation, the Coulomb hole term is reduced 
to 

ScoH(r, rO ^ ^S{t - t'){v{t, t') - W{t, r> = 0)} . 

(35) 

This can be interpreted as the change in the potential as- 
sociated with redistribution of the electron density, which 
is induced by an added electron at the position r. The fac- 
tor 1/2 comes from the adiabatic growth of the screened 
interaction. 

Once the self-energy is computed, the Green's func- 
tion is obtained by solving the Dyson equation (26). Its 
spectral function, ACk^uj), 

A{k,oj) = - V |(^/'k„|ImG(w)|Vk„)| , (36) 

TT ^ — ^ 

n 

is the one electron excitation spectrum for electron ad- 
dition or removal. This is the quantity to be compared 
to (inverse) photoemission measurements. The quasipar- 
ticle energy, which is the peak position of the spectral 
function, is given by 

^k^ = ^kn + Re(V^kn|S(^^^^) - ^xclV^kn) • (37) 

It is often assumed that the self-energy is weakly fre- 
quency dependent. Then the self-energy is safely ex- 
panded around the Kohn-Sham eigenvalue. Taking up 
to the linear order, eq.(37) is reduced to 

^Sf = ^kn + ^knKe(V^kn|S(ekn) - ^xc|V^kn) , (38) 

where the renormalization factor Z is defined by 

^kn = (1 - aReSkn/Sc^|u;=ekJ"' • (39) 

24 LDA + U method 

For the Mott insulator, the LDA and GGA often give 
qualitatively wrong answer. Because of the strong repul- 
sion in short-ranged part of the Coulomb interaction, 
electrons cannot come close each other and the conse- 
quential segregated localization is the origin of the Mott 
insulator. However, LDA treats the distribution of the in- 
teracting electrons as that averaged over the space. Since 
it does not well take into account the electron configura- 
tion avoiding each other, it predicts a metal erroneously. 
On the other hand, if a symmetry breaking such as an- 
tiferromagnetic or orbital order takes place, the electron 
distribution of each spin or orbital component loses its 
uniformity and electrons with different spin or orbital 
avoid each other. If these symmetry breakings are taken 
into account by a mean field, such a segregated local- 
ization can be described. In the Hartree-Fock theory of 
symmetry broken states, though very simple, insulators 
then emerge. 

To fill up and correct the deficiency of LDA, Anisi- 
mov et al. proposed the LDA+U method for the orbitals 
of strongly correlated electrons (typically, d orbital for 
transition metal elements and / orbital for rare earth 
elements), following the spirit of the Hartree-Fock the- 



14 16 ^Yiis method, we start from the Kohn-Sham 
Hamiltonian Hi^j^^x based on DFT and add a correction 
term Ai^ as 

H = Hi^BA + Ai7 (40) 
and then we solve H. Here, the correction is 

AH = JThf - Hoc , (41) 

Hhf = \u X] (42) 

where the density operator rii^j, = <^I,^,f7^«,M,cr at site 
orbital ji and spin a is given by the creation (annihila- 
tion) operator c"*" (c). The Hartree-Fock term i^HF repre- 
sents the short-ranged Coulomb interaction on specified 
electron orbitals fi and v {3d orbital for the transition 
metal elements, for instance and it is denoted as d for 
simplicity) within a unit cell i. (Here we described a sim- 
plified case where the interaction U within the unit cell 
does not depend on the orbitals /i and u. For the moment, 
the exchange interaction is also ignored for simplicity.) 
Since the electron density n^,^ at the unit cell i and or- 
bital fi (including spin degrees of freedom) is determined 
selfconsistently with the solution of Kohn-Sham equation 
after decoupHng, this turns out to be equivalent to the 
level of the Hartree-Fock approximation. 

When we add Hhf^ the Coulomb interaction already 
to some extent counted in the exchange correlation po- 
tential Vxc in LDA is doubly counted. The term i^DC 
is subtracted to remove this double counting. For this 
purpose, 

^DC = \ Y.^''d^{nd^-l) , (43) 

i 

ridi = ^riii, (44) 

is frequently employed. 

The single-particle energy e^^^^ of an electron at the site 
i and orbital u is given from the total energy E = (H) as 
^i,u = d{H) / drii^y. This is rewritten as £^lda + ^(| —^i). 
Then the occupied level at the ith site gives = 1 and 
its energy gets U/2 lower than the LDA energy, while the 
unoccupied level is as large diS U/2 higher than the LDA 
level. This makes the energy difference of U between the 
occupied and unoccupied levels generating a energy gap 
for the electron excitation. This allows a description of 
an insulator. 

The framework can be extended to include the ex- 
change interaction, where eqs. (42) and (44) are modified 
to include the exchange interaction as 

H^^ = l E [(m,m"ic^imV"')<^^'^,";".^"' 

- - {n,n"\U\^"'^'))nl,,-nl^>>.,--- 

(45) 

and the double counting term 

Hj^c = ^ ]^(Undi{ndi - 1) 
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- \AnlMl-^)+niM,,-l)), (46) 

where spin a is separated from the orbital indices ji and 
J specifies the exchange parameter. We introduced the 
notation n^^^, = 

Several different choices of the double counting correc- 
tion Hbc such as fully localized limit and that around 
mean field have been proposed. For details readers are 
referred to Refs.16 and 57. 

In the LDA+U method, one has to determine the value 
of U from the first principles point of view. For this 
purpose, constrained LDA (c-LDA) method is often em- 
ployed. ^^'^^^^ Since U is the quadratic coefficient of the 
n expansion in the Hamiltonian, U is calculated from 
d^E/dii?' by changing n slightly.^^ To change the elec- 
tron density on a specified orbital by keeping the density 
on other orbitals to calculate U on that specified orbital, 
one has to switch off the transfer of electrons between the 
specified and the other orbitals. Then it needs to keep off 
the hybridization contained in the original Kohn-Sham 
equation by hand. This disturbs the original electronic 
structure and introduces an ambiguity in the way of cut- 
ting the transfer. However, if it is done carefully, a rea- 
sonable value of the interaction is obtained. 

In fact, the LDA+/7 method has been applied widely 
to transition metal oxides and has shown reasonable cor- 
rection to the underestimate of the gap known as the de- 
ficiency of LDA.i'^-i^'^2(Q ^ problem with the LDA+/7 
method is that the result may depend on the choice of 
the basis function. In the principle of quantum mechan- 
ics, the ground state should not depend on the choice 
of the basis function. However, since LDA+/7 normally 
takes into account only the short-ranged part of the re- 
pulsion, the final result may depend on the choice of the 
basis function. In addition, the Hartree-Fock approxima- 
tion overestimates the "order parameter" n^,^ as usual as 
one of the mean-field approximation, and ignores quan- 
tum (dynamical) fiuctuations to reduce the localization. 
Then the band gap is usually overestimated. On the con- 
trary, the correlation effect is ignored when the symmetry 
breaking (or localization) is absent. Another drawback is 
the ignorance of spatial fluctuations. A part of the lim- 
itations, namely the missing dynamical fluctuation has 
been examined in a comparison of the static mean fleld 
theory with the dynamical mean fleld theory described 
in the later section. 

3. Downfolding 

3.1 General Framework 

A way to derive an effective low-energy model from 
flrst principles can be performed by calculating the renor- 
malization effect to the low-energy degrees of freedom 
caused by the elimination of the high-energy one, based 
on the idea of the Wilson renormalization group. This 
is the basis of the downfolding method. Choice of the 
low-energy degrees near the Fermi level retained in the 
effective model is not unique. However, if the renormal- 
ization procedure is adequate, flnally obtained properties 
do not depend on the choice. In many of strongly corre- 
lated materials, there exists a well deflned group of bands 
near the Fermi level and isolated from the bands away 



from the Fermi level as we will see examples in Figs. 28 
and 37. This is not accidental because strongly correlated 
electrons at high density as in solid can be realized only 
when the mutual screening of electrons is poor, which 
is ideally seen for isolated small number of bands. This 
isolated group is the natural target of the low-energy de- 
grees of freedom. In transition metal oxides, this group 
consists of the bands whose main component is or- 
bitals at the transition metal atoms. If the crystal fleld 
splitting is large as in many cases with the perovskite 
structure, the low-energy degrees of freedom may fur- 
ther be restricted either to t2g or to eg orbitals only. In 
the organic conductors, the group of HOMO (highest oc- 
cupied molecular orbital) and LUMO (lowest unoccupied 
molecular orbital) are frequently isolated from others as 
we see later. 

In general, a complete set of basis can be constructed 
from the basis functions of a Hamiltonian obtained by ig- 
noring the electron-electron interaction. In this basis, the 
second-quantized total Hamiltonian of electrons equiva- 
lent to the electronic part of eq.(2) is written as 

(47) 

The flrst term represents the kinetic energy including the 
one-body level (chemical potential) term, while the sec- 
ond term is the Coulomb interaction in the present basis 
representation with electronic internal degrees of free- 
dom such as ji. Alternatively one can employ the basis 
of LDA eigenfunctions for Kohn-Sham equation, where 
fi specifles a LDA band and momentum. 

The partition function of this whole electronic system 
is described by Z = Trexp[6'] with S being the action 
given by 

S[c\c] = j Cdr, (48) 

^ = J drc'^drC^H[c'^,c], (49) 

where C is the Lagrangian and x denotes (r, r, a) with the 
spatial coordinate r, the spin a and the imaginary time 
r. Here we have suppressed electronic internal degrees of 
freedom such as the band and momentum index in the 
notation. 

In strongly correlated electron systems, the whole elec- 
tronic Hamiltonians are in many cases rather well sep- 
arated into two parts: One represents electrons in rela- 
tively well isolated bands near the Fermi level and the 
other represents the band degrees of freedom far from 
the Fermi level. Then the Hamiltonian is rewritten in 
the form 

H = Hl^Hh^Hhl- (50) 

In the index, the part representing the bands far from the 
Fermi level (high-energy part) is denoted by the suflix 
while the " "thin-skin" part close to the Fermi level (low- 
energy part) is expressed by the suffix L. The coupling 
between the H and L parts are given by Hhl- Then the 
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trace summation for the partition function is formally 
decomposed to the high- and low-energy parts as 

Z = TTLTiHeMS]- (51) 

After the partial trace over the high-energy part, the 
high-energy degrees of freedom is eliminated leading to 
the action for the low-energy degrees of freedom only as 

SlIcI^cl] = logTr^exp[^]. (52) 

Now the low-energy effective action Sl represented only 
by the L degrees of freedom and the resultant Lagrangian 
Cl = dSL/dr has been derived. When we formally ex- 
pand 

HlIA^cl] = ^l- J c^drcdv (53) 

in terms of the creation and annihilation operators, it 
contains the r dependence in general because of the re- 
tardation effect and the operator Hl is not equal to Hl 
because the partial trace summation over the H degrees 
of freedom renormalizes Hl. The partial trace summa- 
tion may, for instance, be performed by the perturbative 
treatment of the coupling Hhl as we detail later. 

If the dependence on r (namely, the retardation ef- 
fect) can be ignored, Hl may be regarded as an effective 
Hamiltonian. Then the effective Hamiltonian 1-L = Hl in 
general has the form 

n = ^hLK{lJ^,i^)c^^c^, 

+ ^ hLv{lJ^,l~^\T^,y')clc\,CyCy' (54) 

When terms of higher order than this expression (beyond 
the fourth order in the creation and annihilation oper- 
ators) are small, this Hamiltonian is closed within the 
single-particle part and the two-body interactions. If the 
r dependence is appreciable, one has to treat it with the 
action Sl or the Lagrangian. 

To understand the physical applicability of this down- 
folding, let us start from more physically transparent pic- 
ture. Even the electrons belonging to the low-energy part 
originally have the kinetic and interaction energies as in 
the Hamiltonian (47). These energies are, however, sub- 
ject to the renormalization originating from the inter- 
action with the electrons in the high-energy part. This 
effect appears as the "dressing" of the kinetic and inter- 
action energies. Actually, the interaction with the high- 
energy electrons (or holes) Hhl reduces the bare in- 
teraction between low-energy electrons (holes), because 
polarizations of the high-energy electrons (holes) screen 
the interaction between low-energy electrons (holes). In 
addition, for instance, the effective mass is usually en- 
hanced because of the dressing by the high-energy elec- 
trons (holes). This effect is in more general expressed 
by the self-energy effect for the kinetic energy part. The 
real frequency dependence of the screened Coulomb in- 
teraction hLv{^) is obtained from the Fourier trans- 
form and analytic continuation of hLv{^)' In the low 
frequency range that satisfies (1) hLv{^) ^ hLv{^ = ^) 
where the frequency dependence can be ignored, and 



(2) the self-energy approximated as E ~ Rel](cj = 
0) + ujdKeTi/ duj\u;^o^ the retardation effect can be ig- 
nored and the description by the Hamiltonian becomes 
adequate. This corresponds to the case where the low- 
energy electrons can be adiabatically treated under the 
high-energy electrons moving fast. The renormalization 
effect can be ascribed to the screening of the interaction 
and the mass enhancement in the band dispersion given 
by the factor 1/(1 — dReT^ / duj\^=o) multiplying the bare 
dispersion e{k). In examples of the transition metal com- 
pounds, the 3d bands of the transition metal atom are 
relatively well isolated near the Fermi level from others 
as we already mentioned. This makes the description by 
a Hamiltonian appropriate after ignoring the frequency 
dependence of the screening and the mass enhancement 
within the range of the 3d bandwidth. 

3.2 Wannier functions 

To derive the low-energy effective model, one first 
needs to define and specify the low-energy Hilbert space. 
In other words, the first step of the downfolding is to 
construct a set of localized orbitals that span the Hilbert 
space of the low-energy electronic states. In the case of 
SrVOs, for example, three narrow states cross the Fermi 
level (Fig. 5). One may then wish to pick up three lo- 
calized orbitals and construct three-orbital Hamiltonian 
that reproduces the (red) lines crossing the Fermi level 
in the figure. There are several ways for obtaining the 
localized orbitals. Here, we focus on maximally local- 
ized Wannier functions (MLWF) developed by Marzari, 
Souza, and Vanderbilt^^' based on the minimization of 
the quadratic extent of the orbitals. An alternative ap- 
proach is to use the Wannier orbitals of Andersen. The 
former is more general because it does not depend on any 
particular band-structure calculation method. Compari- 
son between the two Wannier functions for some selected 
materials can be found in Ref.66. Although the Wannier 
basis can be chosen arbitrarily in principle and the fi- 
nal results of calculated physical quantities should not 
depend on the choice, it is better to find maximally lo- 
calized orbitals to make the range of transfers and inter- 
actions in the effective lattice models as short as possible. 

Let {V^nk} be the eigenfunctions of the low-energy 
states. Naively, the Wannier function is defined by 

Vnnir) = J e-*-^^/'„k(r)d3k . (55) 

This Wannier function is, however, ill-defined, because 
it depends on the choice of the phase factor at each k 
point. Moreover, at the band-crossing points it is not 
clear which state should be taken. The MLWF utilizes 
this degrees of freedom. The MLWF with band index n 
at cell R is defined by 

^nn{r) = J^J e-^'^-^^il^(r)d3k . (56) 

Here t/j^^ is not the eigenfunction of the Hamiltonian 
(e.g. Kohn-Sham wavefunction) , but it is a linear combi- 
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nation of the eigenfunctions as 

^il^(r) = E^™»(k)^™k(r). (57) 

m 

The coefficients ^/^n(k)'s are numericahy determined 
such that the spread 

^ = ^[{^no\r'^\^no) - {^no\A^nof] , (58) 
n 

is minimized. In contrast with '0nk(r), the gauge of 
'^ik^(^) fixed, and it is a smooth function of k. By 
representing the Hamiltonian in the MLWF basis, 

Hmn{^) = {^mo\H\(PnIl) , (59) 

the on-site energy levels are obtained from m = n, R = 
component. Other matrix elements give the transfer 
integrals. 




Fig. 6. (a) Dyson equation for screened interaction. Thin solid 
lines with arrows represent electron propagators (bare Green's 
functions) Go and thin dashed lines represent bare Coulomb 
interaction V. Bold and bold dashed lines represent their cor- 
responding renormalized Green's function (see Fig. 7) and the 
screened interaction W, respectively. Shaded triangle represents 
the three-point vertex illustrated in (b). (b) Three-point vertex. 




Fig. 5. (Color online) Electronic structure of SrVOs (left). The 
three states crossing the Fermi level (in red) can be treated as the 
low-energy part. The corresponding maximally localized Wan- 
nier function is shown in the right panel. The Wannier function 
having the t2g character is localized around the V atom, with a 
tail at the O sites. 



Fig. 7. Dyson equation for renormalized Green's function. 
Shaded circle represents self-energy illustrated in Fig. 9. Nota- 
tions are the same as Fig. 6. 



^RPA=— ^ — O — 

Fig. 8. RPA diagram of screened interaction, which is obtained 
from W in Fig. 6(a) by ignoring the vertex correction F in Fig. 6. 
This is equivalent to the upper panel in Fig. 4. Notations are the 
same as Fig. 6. 



3.3 Screened Interaction 

Now we discuss how to obtain the renormalization ef- 
fects on the low-energy electrons near the Fermi level 
more concretely. After the partial trace and the elimi- 
nation of the high energy degrees of freedom given in 
eq.(52), the renormalization can be calculated perturba- 
tively in terms of the interaction between the low- and 
high-energy electrons. 

In more general, the Dyson equation for the screened 
interaction is expressed in the diagram illustrated in 
Fig. 6 (a), where the three-point vertex is given in Fig. 6(b) 
and the diagram representation of the dressed Green's 
function is illustrated in Fig. 7. If the vertex correction 
is small as we will discuss later, the screened interaction 
is reduced to the form of the standard RPA shown in 
Fig. 8. In the conventional RPA, we do not distinguish L 
and H degrees of freedom and the polarization (bubble 
in Fig. 6(a) or Fig. 8) contains contribution both from H 
and L electrons. 

Now we extend this conventional RPA to the down- 
folding procedure. In this extension, only the polariza- 
tion containing the high-energy degrees of freedom may 
contribute in the screening and self-energy processes, be- 



cause the polarization purely originating from the low- 
energy degrees has to be removed to keep it dynamical 
in the effective models. Such a RPA with restriction of 
the polarization channel is called the constrained RPA 
(cRPA).^^'^^ In this subsection we sketch how the inter- 
action energy is renormalized by the screening in cRPA. 
In the next subsection we figure out how the kinetic 
energy (namely, dispersion) is renormalized as the self- 
energy effect. 

The bare Green's function is given by eq.(30) while its 
low-energy part Gl restricts the summation within the 
target low-energy bands as 

(60) 

From eq.(60) and below in this and the next subsections, 
the Green's functions are all the bare one but written by 
abbreviating the suffix for a simple notation. 

The total polarization P = Pq in eq.(31) can be di- 
vided into: P = Pl + Ph^ in which Pl includes only Gl 
(i.e limiting the summations in (31) to i, j G {V^l}), and 
Ph be the rest of the polarization. 

Within cRPA, Hlv is reduced to Wh expressed by the 
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bare Coulomb interaction V as 

Wnik) = eJ,\V,k)V{k) , (61) 

ea{w,k) = 1-Pa{k)w{k) (62) 

in the momentum space representation obtained from the 
Fourier transform of r — in eqs.(60), (30) and (31) into 
k, where a represents either L^H or T. Ph is the polar- 
ization function that has contributions from the high- 
energy electrons and k = {k^uo) represents the both de- 
grees of wavenumber and frequency. The RPA form of 
the screened interaction in general is illustrated in Fig. 8 
obtained from the exact form of the screened interaction 
obtained from the Dyson equation in Fig. 6 (a) by ignor- 
ing the three-point vertex correction shown in Fig. 6(b). 
Here, in Fig. 8, cRPA has to contain propagators includ- 
ing some of downfolded electrons, namely "i^" electrons. 
In the lowest order, Ph is given by using the Green's 
function for the low-energy electrons Gl and the total 
Green's function G (or low-energy polarization Pl and 
the total polarization Pt ) obtained from LDA as 

PH(k) = Prik) - PL{k) , (63) 
Prik) = - J dk'G(k')G(k + k'), (64) 

PL{k) = - J dk'GL{k')GL{k + k'), (65) 

where we have divided the whole Green's function G into 
the contribution from the low-energy band, Gl and the 
rest Gh as 

Gik) = GL(k)^GH(k). (66) 

There exists a remarkable identity for the fully 
screened Coulomb interaction W as^^ 

W{k) = e^\V,k)V{k) = el\WH,k)WH{k) . (67) 

This identity proves that the fully screened interaction 
W obtained from the full RPA by the whole polarization 
P is the same as the Coulomb interaction obtained from 
RPA as if one takes Wh were the bare Coulomb interac- 
tion and the low-energy polarization Pl were the full po- 
larization. It assures that one can regard Wh as the effec- 
tive interaction in the effective low-energy model. In this 
cRPA, frequency dependence of Wh has been calculated 
for several transition metal and transition metal com- 
pounds, which revealed that the frequency dependence 
can be ignored within the order of the width of d elec- 
tron bands (typically several eV) (see an example for Ni 
in Fig. 12). This means that the Hamiltonian description 
with the effective interaction U{r) = \\m^^QWH{r^uo) 
becomes appropriate.^^ 

3.4 Self- energy correction 

The renormalization of the kinetic energy is given by 
the self-energy within the level of cRPA. Even when the 
frequency dependence of Wnik) is small in the energy 
range of the low-energy bands, the frequency dependence 
is still large beyond the energy scale of the screening 
channel in cRPA, because Wnik) eventually should come 
back to V{k) in the large uo limit, where the screening 




Fig. 9. Self-energy diagrams expanded in terms of the interaction 
in general. Notations are the same as Fig. 6 and shaded circles 
represent self-energy. Diagrams in the first line in the right hand 
side represent the Hartree-type terms and the second line is for 
the Fock-type exchange contribution. 

does not occur. In the ordinary GW approximation the 
self-energy correction is given by 

T.{k) = j dk'G{k')W{k^k'). (68) 

This is the lowest order term of the self-energy expansion 
represented by the diagrams in Fig. 9, where if G and W 
are replaced with their bare forms. Go and V ^ respec- 
tively, it is reduced to the sum of two diagrams without 
shaded circles in the right hand side in Fig. 9. 

The RPA-level self-energy can be formally rewritten 

as 

E(/c) = j dk'[GL(k')WH(k^k')^GL(k') 

x{W(k^k') - Wnik ^ k')) ^ GH(k)W{k ^ k')]{m) 

In the constrained GW approximation (or equivalently 
cRPA), the first term GL{k')WH{k + k') is excluded be- 
cause this part should be kept dynamical in the effective 
low-energy model. Therefore, the renormalization of the 
LDA band dispersion in the downfolding is given from 
the self-energy correction as 

AS(/c) = j dk'[GL{k'){W{k + k') - WH{k + k')) 

+ GH{k)W{k^k')]-Vxc- (70) 

The last (third) term Vxc is introduced to exclude the 
double counting of the interaction already contained in 
the LDA level. In general the contribution from GhW 
in the second term is small as compared to the first term 
GLiW — Wh)^ because of |Gl| > \Gh\ symbolically and 
from the reason similar to the smallness of the vertex 
correction explained below. 

The band dispersion uj = eo{k) of Hlk obtained from 
LDA is now renormalized by the self-energy AE(/c) ~ 
Eo + Eio; as a; = e%k) = {eo{k) + So)/(l - Si) in the 
low-energy limit. This gives the momentum dependent 
flattening of the dispersions and mass enhancement. The 
imaginary part of the self-energy is normally small which 
validates the hamiltonian description. 

The self-energy correction of the LDA band dispersion 
for the low-energy model has not been extensively stud- 
ied so far and in many cases has been ignored, partly 
because this correction is in general not large. Typically 
the correction reduces the bandwidth of the 3d bands of 
the transition metal oxides as large as 10-20%.^^ 

Another issue of the self-energy effect is related to 
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the double counting already discussed in the LDA+U 
method (§2.4). The interaction effect is already counted 
in the mean-field level in the LDA. It contains the 
Hartree contribution as well as the exchange correlation. 
Therefore, when we consider the self-energy effect, the 
corresponding part of the interaction effect in LDA has 
to be removed. Usually the exchange contribution within 
LDA is expected to be small while the Hartree contribu- 
tion becomes important for the multi-orbital systems, 
because it induces large shifts of the relative level of or- 
bitals. One simple way of subtracting this Hartree contri- 
bution is to adjust the level of orbit als after the Hartree 
aproximation of the effective low-energy model so as to 
have the same level with the LDA result. 

3.5 Low- energy Hamiltonian 

Low-energy effective Hamiltonian in a restricted 
Hilbert space is thus derived by utilizing the hierarchi- 
cal structure of the electronic structure in energy space 
especially applicable to strongly correlated electron sys- 
tems. 

After the Fourier transform and the analytic contin- 
uation from the imaginary time to the real frequency, 
if the frequency dependence of the screened interaction 
Wh{(^) is small and the self-energy can be well approx- 
imated by E(a;) ^ E(a; = 0) -\- ujdReT^/duj\^^o) in the 
energy range of the target bandwidth, the effective low- 
energy model is expressed by a Hamiltonian form of the 
extended Hubbard model, 

H = Hk^Hu, (71) 

Hk = ^ C^nJ^Rn.R'n'CR'n' , (72) 
Rn^R'n' 

Hu = - ^ C^Rn^Rn'Unn'R.rnm'R'C^R'rn^R'rniJ'^) 
R,nn' ,mm' 

after the cRPA procedure. Here, n, n', m, m' denote both 
of spin and orbital degrees of freedom and R is the spatial 
coordinate. Now we have arrived at the effective low- 
energy Hamiltonian by the downfolding to be solved by 
low-energy solvers discussed in §4. 

3.6 Vertex correction 

In the standard downfolding scheme, the cRPA is an 
efficient way to derive the renormalized kinetic and inter- 
action energies in the effective models. Usually the con- 
ventional RPA is justified for the case \PT{k)V{k)\ < 1 
in the denominator of eq.(61) with eq.(62) in the per- 
turbative sense. However, in the present case, the bare 
Coulomb interaction is typically as much as several tens 
of eV, while the polarization of the metallic target band 
is scaled by Pt ~ A£^~^, where AE is the typical en- 
ergy scale of the target bandwidth, because Pt(^) ^ 
J dk'G{k)G{k + k') is basically scaled by the energy de- 
nominator of the Green's function for the downfolded 
bands as eq.(30), and G is given by eq.(60). Then AE is 
typically the inverse of eV. Therefore, \PT{k)V {k)\ ^ 1 
typically holds and clearly violates the above requirement 
meaning that the simple full RPA cannot be used. 

Even in the case of cRPA, the polarization contain- 



(a) rio= J>- (b) Ti - J>— 

Fig. 10. Three-point vertex diagrams, (a) Lowest order vertex (b) 
Dressed vertex represented with the renormahzed interaction and 
Green's function. Notations are the same as Fig. 6. 

ing the high-energy bands {H bands) is scaled by P ~ 
A£^~^, where AE is now the typical energy scale of the 
downfolded band {H bands) measured from the Fermi 
level, because Ph is scaled by the energy denominator of 
Gh for the downfolded bands. Then AE is typically the 
inverse of several eV. Therefore, \PH(k)V {k)\ ^ 1 still 
holds. 

Nevertheless, even in this region, cRPA can be a good 
and convergent approximation as we describe below: We 
can show that cRPA becomes accurate and convergent 
when the vertex correction is small (|r — 1| <C 1), namely, 
the difference between the three-point vertex diagram 
F shown in Fig. 6(b) and its lowest order term Fq = 1 
(the first term in the right hand side) is small, because 
the dressed interaction with the full vertex correction 
W shown in Fig. 6 (a) then becomes nearly the same as 
the RPA diagram in Fig. 8. Here in cRPA, at least some 
of the propagators have to contain the electrons in the 
downfolded band ("i^" electrons), because the diagram 
containing only the "L" electrons should be excluded in 
cRPA. 

Let us consider the lowest order term of the vertex 
correction in Fig. 10(a) Again, this lowest-order term is 
basically proportional to PhV ~ V/ AE^ which is much 
larger than unity. On the other hand, when we con- 
sider the dressed vertex correction shown in Fig. 10(b), 
the renormalized interaction vertex and the renormal- 
ized propagators should be employed in a self-consistent 
fashion. Then the renormalized interactions between two 
H electrons and/or between a L electron and a H elec- 
tron appear, together with propagators consisting of two 
H electrons G ^ ipni^H and/or of H and L electrons 
G ~ i^ni^L' Therefore, the vertex correction is expanded 
in terms of W^/AE, where is the interaction ei- 
ther between two H electrons or between a H electron 
and a L electron, screened by H electrons. Then 
is in general even smaller than Wh in eq.(61) because 
of more efficient screening by the electrons on the same 
band for the former and because of the interaction at 
distant L and H Wannier orbit als in the latter. There- 
fore, W^fj/AE can in general be a small parameter. In 
addition, if G ^ i^ni^L is contained, it multiples another 
small factor as a small matrix element in the numerator 
of the vertex correction, because the overlap of the wave- 
function on different bands and iIjh is small. These 
are in general justified if the coupling between the target 
and the downfolded bands are weak, in either sense of 
the hybridization ( or wavefunction overlap) and/or the 
interband interaction. 

Even more important reason for the irrelevance of the 
vertex correction in cRPA is that in the counting of W'^ , 
the full screening channel by the L electrons should be 
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included, though this channel is excluded in the counting 
of Wh' The gapless particle- hole excitations of the L 
electrons very efficiently screen the Coulomb interaction 
as is known in the full RPA results shown in Fig. 12.^^ 
In short, the vertex correction is irrelevant provided that 
W^/AE is small and this criterion is not perfectly but 
rather well satisfied in real strongly correlated electron 
materials. In typical transition metal oxides, may be 
less than 1 eV while AE is several eV. 

Although it may be small, but the vertex correction 
arising from the W^^/AE expansion quantitatively en- 
hances the screening. On the other hand, the polariza- 
tion from the transition between L and H electrons may 
be quantitatively reduced when the L electrons are under 
strong correlation, because the energy for adding an L 
electron and removing an L electron may split as is typ- 
ical in the Hubbard band splitting of the L band to the 
upper and lower Hubbard bands. This splitting will re- 
duce the screening and at least partially cancel the above 
vertex correction. Therefore, the vertex correction may 
be even smaller. 

The self-energy coming from the H electrons discussed 
in §3.4 becomes also small from the same reason in the 
same situation. These small vertex correction and small 
self-energy are the reason why the cRPA offers a good 
approximation. The irrelevance of the vertex and self- 
energy correction is understood more intuitively; In the 
insulator or semiconductor where the density of states is 
zero around the Fermi level, the vertex and self-energy 
correction become small when the interaction is smaller 
than the gap. The renormalization by cRPA is similar, 
where the gapless particle-hole excitations are excluded. 
Though the correction is small, there must exist some 
finite correction from the vertex part, while its quantita- 
tive estimate has not been fully examined so far and is 
left for future studies. 

3. 7 Disentanglement 

Although the cRPA method offers a general and accu- 
rate scheme for the downfolding, its applications to real 
systems still have a serious technical problem. The prob- 
lem arises when the narrow band is entangled with other 
bands, i.e., if it is not completely isolated from the rest 
of the bands, which is the case in many materials. Even 
in simple materials such as the 3d transition metals, the 
3d bands mix with the 4s and 4p bands. Similarly, the 4/ 
bands of the 4/ metals hybridize with more extended s 
and p bands. For such cases, it is not clear anymore which 
part of the polarization should be eliminated when calcu- 
lating the screened interaction using the cRPA method. 

In this subsection, we take the notation of "d" and 
"r" symboHcally, instead of L and H to specify the 
low-energy target bands and the high-energy downfolded 
bands, respectively employed in the previous sections, to 
deliver a more concrete image of the d bands and the 
rest bands in transition metals and transition metal com- 
pounds. Since it does not mean the loss of generality, one 
may interchange the notations each other. 

If the d subspace forms an isolated set of bands, as for 
example in the case of the t2g bands in SrVOa as we saw 
in Fig. 5, the cRPA method can be straightforwardly ap- 



plied. However, in practical applications, the d subspace 
may not always be well identified. An example is 3d tran- 
sition metal series such as Ni shown in Fig. 11 (a), where 
the 3d bands are entangled with the 4s and 4p bands. 

To treat this complexity, a prescription for entangled 
bands was proposed recently.^^ The essential point is 
that one has to strictly keep the orthogonality between 
the low-energy subspace contained in the model and the 
complementary high-energy subspace to each other. The 
orthogonalization by the projection technique enables a 
proper disentanglement of the bands. Although physi- 
cal properties may not sensitively depend, we still have 
a freedom that the d space somehow depends on the 
choice of the energy window when one constructs the 
Wannier functions. However, once the disentangled band 
structure is obtained, the constraint RPA method can be 
used to determine the partially screened Coulomb inter- 
action uniquely. Numerical tests for 3d metals show that 
the method is stable and yields reasonable results. The 
method is applicable to any system, and applications to 
more complicated systems. We review this procedure in 
more details here. 

We first construct a set of localized Wannier orbitals 
from a given set of bands defined within a certain en- 
ergy window. These Wannier orbitals may be generated 
by the post-processing procedure of Souza, Marzari and 
Vanderbilt^^' or other methods, such as the preprocess- 
ing scheme proposed by Andersen et at. within the Nth- 
order muffin-tin orbital (NMTO) method. We then fix 
this set of Wannier orbitals as the generator of the d 
subspace and use them as a basis for diagonalizing the 
one-particle Hamiltonian, which is usually the Kohn- 
Sham Hamiltonian in LDA or in generalized gradient 
approximation (GGA). The obtained set of bands, defin- 
ing the d subspace, may be slightly different from the 
original bands defined within the chosen energy window, 
because hybridization effects between the d and r spaces 
are switched off. It is important to confirm that the dis- 
persions near the Fermi level well reproduce the original 
Kohn-Sham bands. 

The wavefunctions are thus projected to the d space 

by 

1^.) = V\^^) , (74) 
where the projection operator V is defined as 

We define the r subspace by 

|0,) = (1-P)|^,) (76) 

which is orthogonal to the d subspace constructed from 
the Wannier orbitals. In practice it is convenient to 
ortho normalize and prepare N — Nd basis func- 

tions. By diagonalizing the Hamiltonian in this subspace 
a new set of wavefunctions {0^} and eigenvalues {e^} 
{i = 1, • • • — Nd) are obtained. Namely, the Kohn- 
Sham Hamiltonian becomes block diagonal in the d space 
and r space separately, and the hybridization effects be- 
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f d space 
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r space 



(77) 



As a consequence of orthogonalizing {cpi} and {V^j}, the 
set of r bands {e^} are completely disentangled from 
those of the d space {sj}^ and they are slightly differ- 
ent from the original band structure {si}. As we will see 
later, however, the numerical tests show that the disen- 
tangled band structure is close to the original one. 

From the "d" bands, we calculate the "d" polarization 
Pd as 



occ unocc 



P,(r,r';^) = ^^ 



V'*(r)^/;,(r)V';(r')V'i(r') 



IT] 



IT] 



(78) 



where {V^i}, {si} {i = 1,---A^c^) are the wavefunctions 
and eigenvalues obtained from diagonalizing the one- 
particle Hamiltonian in the Wannier basis. 

The effective screened interaction for the low-energy 
model is calculated according to eq.(62) with Pr = 
P — Pd^ where P is the full polarization calculated for 
the disentangled band structure. It is important to re- 
alize that the screening processes from the polarization 
Pr include the Coulomb interaction between the d space 
and the r space, in calculating the screened interaction 
of d bands (so called U terms in eq.(73)), although the 
d-r hybridization is cut off in the construction of the 
wavefunctions and eigenvalues. 

We also note that starting from the complete orthog- 
onality between the d and r bands in Eq.(77) is crucially 
important to assure stable cRPA calculations. In fact, if 
the orthogonality is not perfect, the resultant frequency- 
dependent screened Coulomb interaction could have un- 
physically negative values in some frequency region. 
For example, if Pr = P — Pd would be calculated from 
Pd obtained from the above procedure while the total P 
is calculated from the original LDA band, the resultant 
Pr gives unstable behavior of the screened interaction /7, 
because such a Pr is not compatible with mutually or- 
thogonal subspace of r and d and a small nonorthogonal- 
ity between d and r subspaces implicitly assumed in this 
Pr would yield a singular behavior at low energies. One 
has to use the total P obtained after the disentanglement 
procedure to assure the orthogonality of the subspaces. 

This disentanglement procedure has been tested to 
work in 3d transition metals. For technical details see 
the literature. ^^^^ Figure 11(a) shows the Kohn-Sham 
band structure of nickel. There are five orbitals having 
strong 3d character at [-5 eV:l eV], crossed by a disper- 
sive state which is mainly of 4s character. Using the pre- 
scription for the maximally localized Wannier function, 
and with the energy window of [-7 eV:3 eV], interpolated 
"d" bands are obtained. The subsequent orthogonaliza- 
tion procedure gives the orthocomplementary "r" bands. 
Comparing Fig. 11(b) with (a) we can see that there is 
no anti-crossing between the d bands and the r bands in 
(b) contrary to (a). Aside from this difference, the two 
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Fig. 11. (Color) (a) Kohn-Sham band structure of nickel in LDA. 
(b) Disentangled band structure with d-r hybridization switched 
off. The red lines show the d states obtained by the maximally 
localized Wannier scheme, while the blue lines are disentangled 
r states. Energy is measured from the Fermi level. 



band structures are very similar. 

The effective screened interaction for the effective 3d 
model is calculated by the constrained RPA, namely, by 
eq.(61). The results are shown in Fig. 12. First, we observe 
that the frequency dependence of U is small within the 
3d bandwidth (~ 5 eV), which justifies the treatment by 
an effective model Hamiltonian within this energy range. 
Second, as is expected, the partially screened onsite inter- 
action U from eq.(61) (~ 4eV) is significantly larger than 
the fully screened RPA result W from eq.(67) {^1 — 2 
eV). At low frequencies (namely, for < 5 eV). This im- 
plies proper elimination of d-d screening processes has 
a large effect. This was also applied to a series of other 
3d transition metals and was found to give stable and 
reasonable results. 

In the present formulation, small off-diagonal matrix 
elements of the Kohn-Sham Hamiltonian between the d 
wavefunction \ipi) and the r space are ignored. How- 
ever, if the energy of the d-r hybridization point in the 
band dispersion is smaller than the energy scale of inter- 
est, one has to retain all of these hybridizing bands in the 
effective model, because the hybridization effect changes 
the band dispersion and the wavefunction significantly 
in the vicinity of the anti-crossing points. In many cor- 
related materials with d or f electrons, the energy scale 
of interest determining material properties is typically of 
the order of 100 K or lower, which is smaller than the 
typical energy crossing points. Therefore, the low energy 
models constructed only from the dor f Wannier orbitals 
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Fig. 12. (Color online) Effective onsite Coulomb interaction U 
calculated by cRPA (filled (red) circles) (eq. (61)) and fully 
screened onsite interaction W by RPA (open (blue) circles) 
(eq.(67))for nickel as a function of frequency by disentanglement 
procedure. The diagonal terms of screened interactions aver- 
aged over the 3d orbitals in the Wannier basis are plotted. 



may give at least a good starting point of understanding 
the low energy physics. 

3.8 Dimensional downfolding 

Remarkable progress in understanding physics of ma- 
terials with low-dimensional anisotropy such as cuprate ''^ 
and iron-based^^ superconductors have stimulated stud- 
ies on electronic models in low dimensions. In particu- 
lar, ID or 2D simplified models are frequently used and 
have greatly contributed in revealing characteristic low- 
dimensional physics with strong-correlation and fluctua- 
tion effects.^ 

However, the ab initio downfolding method reviewed 
in this article are formulated to derive ab initio models 
in 3D space. Thus, we have to solve the derived model 
as a 3D model, but this requires significantly demanding 
computation when we consider correlation effects with 
high accuracy. To construct a low-dimensional model 
tractable by the widely employed theoretical approaches, 
it is crucial to bridge the 3D models to effective ID or 
2D models based on the ab initio derivation. 

Recently, a scheme of downfolding a 3D model to 
lower-dimensional models from first principles has been 
formulated as a dimensional downfolding in real space. 
This supplements the original band downfolding in en- 
ergy space. The formalism eliminates the degrees of free- 
dom for layers (chains) other than the target layer (chain) 
after the interlayer /chain screening taken into account. 
This is useful when low-energy solvers for the effec- 
tive models allow only the 2D models as tractable. The 
scheme is general and particularly works well for quasi- 
low-dimensional systems. The dimensional downfolding 
has another computational advantage that the range of 
the screened effective interaction becomes short-ranged 
when we take into account interlayer (interchain) screen- 
ing for metals. Note that, the original band downfolding 
reviewed in previous sections leaves the effective screened 
interaction long-ranged even for metals because cRPA 
excludes metallic screenings in the target bands. 

The polarization in the target band is decomposed 



into layer-by-layer contributions in the real space. The 
RPA polarizations except for those within the target 
layer/chain (namely the processes in Figs. 13(b) and (c) 
excluding the process in Fig. 13(d)) contribute to the 
interlayer /chain screening, which renormalizes the ef- 
fective interaction between electrons within the target 
layer/chain. This screening deletes the long-range part of 
the interactions for the case of metallic systems and justi- 
fies the short-ranged models as effective ab initio models 
of real materials. 

As an example, it was applied to derive an effective 
2D model for LaFeAsO. It was found that the inter- 
layer screenings reduce onsite Coulomb interactions by 
10-20 % and further remove the long-range part of the 
screened interaction. This formalism justifies a multi- 
band 2D Hubbard model for LaFeAsO from first prin- 
ciples. 




Fig. 13. (Color online) Schematic diagram for effective interac- 
tions between electrons at ri and r2 in the target layer/chain, 
screened by intra- and inter-layer/chain polarizations x*(r,r^). 
(a) shows general second-order diagram for the screened inter- 
action (b) shows the screening by a polarization in the other 
layers/chains, while (c) describes that by an interlayer polar- 
ization between the target and the other layers, (d) shows the 
screening by a polarization within the target layer/chain itself. 
This process should be excluded in the present 2D-cRPA down- 
folding, while (b) and (c) are included. Notations are the same 
as Fig.6.'^^ 



Here, we discuss vertex corrections and self-energy ef- 
fects. The band downfolding (3D-cRPA) becomes ade- 
quate basically when the downfolded band energy in the 
denominator of the propagator is far from the Fermi level 
as we discussed in §3.6. In the case of the dimensional 
downfolding, the 2D-cRPA (or ID-cRPA) becomes ac- 
curate not by a large energy denominator but by small 
numerators of the vertex expansion. If two neighboring 
layers or chains are far apart in space, the overlap of the 
wavefunction between two Wannier orbitals on different 
layers/chains are small. In this case, the Green's function 
(60) with one i/j being on one layer/chain and its part- 
ner being on another layer/chain (symbolically shown in 
Fig. 13(c) for the lowest order) hardly contributes to the 
polarization for the screening channel. Even in this case, 
contribution of the polarization from the propagator with 
both of t/j on the off-target layer/chain is large. With this 
polarization, the 2D-cRPA diagram contains screened in- 
terlayer /chain interaction to connect the target and off- 
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target layer/chain (see Fig. 13(b)). This screened inter- 
layer /chain interaction Winter is again given as a con- 
sequence of efficient screening by the gapless particle- 
hole excitations on the target and oflF-target layers and 
it is normally smaller than discussed in §3.6. Since 
the vertex is roughly given by the expansion in terms 

of Pl Winter ^ Winter/ Wintra, whcrC Wintra IS a typical 

screened intralayer/ chain interaction calculated from the 
full RPA. The ratio of these two interactions turns out to 
be the small parameter in the vertex correction for the 
dimensional downfolding. The self-energy contribution in 
the dimensional downfolding has a similar feature where 
we have Winter as a small parameter in the expansion. 

Here, we note a more fundamental issue. Of course, 
the present dimensional downfolding does not mean that 
the 3D systems can be rigorously mapped to lower- 
dimensional models. For example it is obvious that we 
are unable to treat the 3D ordering process. Nevertheless, 
within the properties and questions that allow the neglect 
of the energy scale of the interlayer / chain effective trans- 
fer and effective interaction, the present downfolded low- 
D models offer the best way of simplifying the problem to 
that for low-dimensional physics in an ab initio way. An 
important point is that the interlayer / interchain interac- 
tion is more efficiently screened in the band downfolding 
procedure and can be much smaller than the original ef- 
fective interaction for 3D systems. 

4. Low-Energy Solver 

The next task of the three-stage scheme is to solve the 
effective models.^ Historically, strongly correlated elec- 
tron systems have been studied extensively by using an 
ad hoc theoretical models. Such theoretical models con- 
tain only small number of bands and restricts the inter- 
action range short as in the cases of the Hubbard model, 
Anderson model and Heisenberg model. Such simplifi- 
cations have made it possible to carry out high-accuracy 
calculations including ordering and quantum fluctuations 
beyond the mean- field level employed in DFT. In fact, 
mechanisms of antiferromagnetic phase stabilized by the 
superexchange interaction, Kondo effect, and Mott tran- 
sitions have been elucidated by using theoretical models. 
In those studies in the long history, the interaction pa- 
rameter U in the Hubbard model has been chosen by 
hand to satisfy physical intuitions and/or comparisons 
with experimental results. However, rapid progress in re- 
search for transition metal compounds including the cop- 
per oxide superconductors and rare earth compounds dis- 
playing heavy-fermion behavior as well as discoveries of 
many correlated and functional electron systems in mate- 
rials research have strongly promoted studies of deriving 
models of real materials based on the first-principles cal- 
culations. This trend is supported by the fact that one 
can not have a clue for understanding mechanisms of in- 
triguing phenomena when an approximate solution of an 
ad hoc model does not reproduce experimental results; 
the disagreement could be ascribed to poor approxima- 
tions in solvers, while it could equally be ascribed to a 
false of the model itself. The present three-stage RMS 
scheme with the downfolding procedure has opened an 
avenue of studies on theoretical models on a firm basis 



of first principles to overcome such uncertainties. 

Several different low-energy solvers have been applied 
to solve effective lattice Fermion models derived by the 
three-stage RMS scheme. The present low-energy solvers 
are roughly classified into two streams. One is the meth- 
ods for solving lattice Fermion models as quantum many- 
body problems. The other is the methods based on dy- 
namical mean-field theory. 

Dynamical mean field theory (DMFT) is formulated 
by ignoring the momentum dependence of the self-energy 
but considering the frequency dependence correctly. 
This becomes a good approximation when the spatial 
dimension increases and is proven to be exact in infinite 
dimensions. In other words, the DMFT becomes exact 
when the coordination number becomes infinite. 

For lattice Fermion solvers, various methods as exact 
diagonalization, auxiliary field Monte Carlo (AFMC),^^ 
path-integral renormalization group (PIRG),^^'^^ den- 
sity matrix renormalization group (DMRG), many- 
variable variational Monte Carlo (mVMC)^^ and 
Gaussian-basis Monte Carlo^^'^^ have been developed, 
which are applicable when the models have the Hamil- 
tonian expression eq.(71). In principle, these approaches 
allow taking account of spatial fluctuations equally with 
dynamical fluctuations in contrast to DMFT, while com- 
putational cost becomes demanding. Functional renor- 
malization group (fRG) has also been developed. ^^^^ 
This method divides the whole Brillouin zone into 
patches and the renormalizations of the coupling con- 
stants in each patches are calculated with renormaliza- 
tion group transformation by approaching the Fermi en- 
ergy window. Since it gives how various coupling con- 
stants to orderings grow and which instability occurs 
first, it is suited in the weak-coupling region. Within 
more biased framework, or within weak-coupling or 
mean- field framework, simple RPA or fluctuation ex- 
change approximation (FLEX) have also been used as 
computationally tractable methods. Below we review 
the DMFT (either combined with GW or without it), 
mVMC and PIRG as typical and extensively examined 
solvers . 

4.1 Dynamical mean- field theory 

Usual mean-field theories approximate effects of in- 
teracting electrons with a static effective field. As a re- 
sult, the many-electron problem is reduced to a single- 
particle problem under the influence of the static mean 
field. In the dynamical mean field theory (DMFT),'^^'^^ 
the many-electron problem is mapped onto an impurity 
problem in the following way. A chosen site in a periodic 
lattice is treated as an impurity and the surrounding sites 
are treated as a bath. The effective Coulomb interaction 
between electrons at the impurity is taken into account 
explicitly, whereas the self-energy arising from the sur- 
rounding bath is taken into account as a dynamic mean 
field. The impurity self-energy Sinip(cj) is energy depen- 
dent but is assumed to be local, i.e., not k dependent. 
Then, the self-energy of the original lattice model is re- 
placed by the impurity self-energy on each site. Thus, 
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the lattice Green's function is given by 

1 



The impurity Green's function Gimp 

G™p(t) = -(Tco(t)c+(0))^^^^ 
is calculated from an effective action, 



(79) 
(80) 



^eff 



JdrJ dr'Y,ctir)go'{r-r')c^Ar') 



(81) 



Here the subscript in eqs.(81) and (80) denotes the im- 
purity site and co(r) is a Grassmann variable. The band 
index is omitted for simplicity. The dynamical mean field 
Qq^ is given by 

gQ^{iujn) = G^^l{iuJn) + Simp(^a;n) , (82) 
where the local Green's function is defined by 

1 



Gloc(^^n) 



/Li - ifo(k) - Simp(iw„) ' 



(83) 



^ — ,L , — / — 'impV 

To understand the meaning of Qq^ consider the following 
expression: 

a(k)-i =G(k)-i+E,^^. (84) 

Since G is the full Green's function defined in eq.(79), Q 
is the lattice Green's function excluding the effect of the 
self-energy at the impurity site. Qo is the projection of 
Q at the impurity site which contains the effects of the 
self-energy from the rest of the sites. Consequently Qo is 
different from the non-interacting Green's function. 

The self-consistency condition requires that the local 
Green's function is equal to the impurity Green's func- 
tion, 

Gloc(^^n) = Gimp(^^n) • (85) 

Thus the self-consistent calculation is carried out in the 
following way.^'^^ 

(1) For a given mean field Qq^^ the effective action 6'eff 
is determined. 

(2) The impurity problem eq.(81) is solved, and the im- 
purity Green's function Gimp and the impurity self- 
energy Eimp = Go^ - G^J^^ are calculated. 

(3) The impurity self-energy is used for the lattice self- 
energy. The lattice Green's function is computed 
from eq.(79). 

(4) The local Green's function is calculated from 
eq.(83). 

(5) If the self-consistency condition eq.(85) is not sat- 
isfied, a new mean field Qq^ is constructed from 
eq.(82), and the self-consistent cycle is continued. 

There are several techniques for calculating the impu- 
rity Green's function eq.(80) from eq.(81), such as it- 
erative perturbation theory (ITP), exact diagonaliza- 
tion, numerical renormalization group and various quan- 
tum Monte Carlo methods. Numerical renormalization 



group offers a high accuracy method at low energies. 
A quantum Monte Carlo method with the Hirsch-Fye 
algorithm^^ is widely used for obtaining finite temper- 
ature properties. In this algorithm, the imaginary time 
in the path integral must be discretized, which necessi- 
tates an additional extrapolation procedure for the con- 
tinuum limit. Recently by utilizing a continuous time al- 
gorithm,^^' an efficient solver has been developed in 
the weak coupling approach^^'^^ as well as in the strong 
coupling approach, where the exponential of the Hamil- 
tonian is expanded^^'^^ in terms of either the interaction 
or the kinetic energies and it is free of the discretiza- 
tion error. A variation called diagrammatic Monte Carlo 
method has also been proposed. The implementation 
of nonlocal correlation effects beyond DMFT has also 
been attempted from other approaches such as the dual 
Fermion method^^'^^ and dynamical vertex approxima- 
tionioO'ioi 

The DMFT successfully describes the correlation- 
driven metal-to-insulator transition in one consistent 
theoretical framework as we see the evolution of the den- 
sity of states in Fig. 14. As U/t increases from a weakly 
correlated metallic regime, the quasiparticle peak gets 
narrower. At the same time, the upper and lower Hub- 
bard bands evolve. As U/t increases further, the quasi- 
particle peak disappears and the system becomes gapful. 




Fig. 14. Density of states A{ou) = -ImG by single-site DMFT at 
T = for half-filled Hubbard model at U/W = 0, 0.2, 0.4, . • • • 
and 1.6 from top to bottom around oo = 0, where W is the 
noninteracting bandwidth. The original noninteracting density 
of states is taken as semicircular. The calculation is done by 
numerical renormalization group. At U/W = 1.6 is an insulator 
to which the coherent peak in metals becomes sharpened until 
U/W = 1.4.102 See also Fig.l5(a) 



To improve the lack of k (momentum) dependence of 
the self-energy in the single-impurity DMFT, there are 
several attempts to include the k dependence in the self- 
energy. Most commonly used methods are the dynami- 
cal cluster approximation^^^ and cellular DMFT,^'^' ^^"^ in 
which the impurity problem is defined not for a single site 
but for a cluster including several sites (or alternatively 
Brillouin zone is divided into patches to allow the mo- 
mentum dependence of the self-energy). Cluster pertur- 
bation theory was developed by Senechal et al. to include 
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intercluster coupling as an RPA type perturbation. 
Potthoff formulated a Baym-Kadanoff-type formalism as 
a functional of the self-energy instead of Green's func- 
tion in case of DMFT/^^ It was applied to a formalism 
for cluster degrees of freedom called variational cluster 
approach, which bridges DMFT and the cluster per- 
turbation theory. 

A typical effect of spatial correlations is seen in the 
cluster extension of DMFT for the Hubbard model. It 
shows a suppression of the coherent peak with a pseudo- 
gap formation in contrast to the sharpening of the coher- 
ent peak at the Fermi level close to the metal-insulator 
transition in the single site DMFT (see Fig. 15). This 
indicates an appreciable role of intersite (spin) correla- 
tions ignored in the single site study depending on the 
lattice structure. Differentiations of electrons in momen- 
tum space make the Mott transition strongly momentum 
dependent with distinction of more correlated and less 
correlated regions. On the square lattice, the pseudogap 
opens first in the "antinodal" region (around (tt, 0) and 
(0,7r) regions) in metals separated from the opening of 
the real Mott gap eventually around the "nodal" region 
at the Mott transition. ^^^^^"^ 



1-siteGDMFT 
- - - U=t 




Fig. 15. (Color online) Density of states of Hubbard model on 
half-filled square lattice obtained from CDMFT.^*-*^ Results from 
the single-site DMFT (a) and 4-site cellular DMFT (b) are com- 
pared for four different choices of U jt. The broadening factor 0.1 
is employed. In contrast to (a), pseudogap is developed in (b) 
near the transition to the insulator. 



A combination of DMFT and LDA, so-called 
LDA+DMFT, was developed to apply DMFT to real ma- 
terials.^^' ^^^^^^ As is the same as the LDA+U method, 
LDA+DMFT is based on a many-body Hamiltonian, 



E 

{imcr} 



LDA 

imA' m' 



imm' a 



where i^^DA ^.j^^ Kohn-Sham Hamiltonian in the LDA. 
Hubbard terms for direct and exchange interactions for 
the correlated orbitals, e.g. d oi f orbitals, are added on 
top of the LDA Hamiltonian. In order to avoid double 
counting of the Coulomb interactions for these orbitals, a 
correction term i^DC is subtracted. The resultant Hamil- 



tonian (86) is solved by the DMFT by assuming that 
the many-body self-energy associated with the Hubbard 
interaction terms can be calculated from a multi-band 
impurity model. As described above, the method can be 
applicable for a wide range of /7, from metallic regime, 
such as Fe and Ni, to strongly correlated insulator, 
such as NiO, including the intermediate regime where 
both the coherent and the incoherent peaks exist (e.g. 
SrVOa). On the other hand, the first-principles determi- 
nation of U and proper treatment of the double counting 
term are challenges as is discussed in this review in detail. 

A combination of the GW approximation and DMFT 
is a new approach to go beyond the LDA+DMFT 
method. They were proposed both in a model context^^^ 
and within the framework of realistic electronic struc- 
ture calculations. The basic idea is simple. The DMFT 
is suitable to treat the onsite correlations, while the 
k dependence in the self-energy and long-range inter- 
action effects are not taken into account. The RPA, 
on which the GW approximation is based, is generally 
good for handling long-range correlations. Hence, in the 
GW-hDMFT approach, the on-site self-energy is taken to 
be the DMFT self-energy, while the off-site self-energy 
is calculated by the GW approximation. Viewed from 
the GW, the on-site GW self-energy is supplemented by 
that of DMFT, correcting the GW treatment of on-site 
correlations. Viewed from the DMFT, the off-site contri- 
butions to the self-energy approximated within the GW 
approximation give a momentum dependent self-energy. 

The above idea can be formulated using the following 
free-energy functional, 

T{G,W) = TrliiG -Tr[{GH^ - G-^)G] - \Tr\iiW 



-Tr[{v-^ - W-^)W] + ^[G, W] , (87) 



where G^^ 



iujn + /i + /2 — Vh is the Hartree Green's 
function with Vh being the Hartree potential. The func- 
tional is an extension of the functional by Luttinger and 
Ward (LW),^^^ and has two variables, the Green's func- 
tion G and the screened Coulomb interaction W. By tak- 
ing functional derivatives of ^ with respect to G and 
the stationary condition 



ST _ ST 
SG~ ' SW 







yields 



G-^ = c^^-s, j: = s^/6g, 

W-^ = V-^-P, P = -26^ /5W. (89) 



Now, the functional ^ is divided into two parts as 



where R denotes a lattice site. The first term is approx- 
imated in the GW approximation as 



<i>Gw[G,W] = ^GWG 



(91) 



while the second term is evaluated from the impurity 
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problem defined by the following action, 

S = I drdr' [-^ct{r)gzl,{T-T')cL>{r') 

x:<(r')cL,(r'):], (92) 

where the double dots denote normal ordering and L 
is an orbital of angular momentum on a given sphere 
where the impurity problem is defined. Then the above 
stationary conditions yield 

P = Pgw{^- ^RR' ) + PiZ^RR' • (93) 
In momentum space, eqs.(93) are expressed as 

(94) 
(95) 

The outline of the self-consistent loop is the following. 

(1) The impurity problem (92) is solved for a given 
Weiss field Qll' and interaction Uais^ and the im- 
purity Green's function Gimp and self-energy Eimp 
are obtained. The two-particle correlation function 

XL^L^L:,L^ = (: c[^{r)cL2{r) :: c[^{r)cLA^') :)s 

(96) 

is also evaluated. 

(2) The polarization function of the impurity problem 
is computed from the interaction Uai3 and the cor- 
relation function eq.(96). 

(3) The full Green's function G(k, iun) and effective in- 
teraction W(q, ii^n) are constructed from eqs. (94) 
and (95). Their local part is defined by 

Gloc(^^n) = ^G{k,iujn) , (97) 
k 

Wiociii^n) = ^Wiq.iun) . (98) 
q 

(4) The Weiss dynamical mean field Q and the interac- 
tion U are updated according to 

= G^ol + Simp , (99) 
= Woe + Pi^P ■ (100) 

This cycle is iterated until self-consistency for Q and U is 
achieved. When self-consistency is reached, Gimp = Gioc 
and Wimp = W\oc are satisfied. Therefore, the second 
term in eq.(94) can be rewritten as 

^S^^(k,r) = - ^ wtJ:f''hr)Gfi^/iT) (101) 

k LiL[ 

This shows that the on-site contribution of the GW self- 
energy is precisely subtracted out, thus avoiding double 
counting. 



Eventually, self-consistency over the local electronic 
density can also be implemented. Using the new density 
from the Green's function at the end of the convergence 
cycle above, the next iteration of the GW calculation can 
be carried out, until the self-consistency with respect to 
the charge density is achieved. 

So far, the full implementation of the above scheme is 
not yet done. Instead, a simplified scheme was applied 
to nickel. In the application, the GW calculation was 
done only for one-shot. Also, the frequency dependence of 
U was neglected, and its static value was used. Improve- 
ment over these simplifications and more applications are 
highly anticipated. 

4-2 Variational Monte Carlo method 

Many-variable variational Monte Carlo method has 
been developed recently as a choice of the low-energy 
)^olver.^^' Historically, the variational wavefunction 
has played important roles in understanding physics 
of interacting fermions. Superconductivity by the BCS 
wavefunction, the roton excitation of ^He by Feynman^^^ 
and Laughlin's wavefunction for the fractional quantum 
Hall state^^^ are typical examples. 

However, it is also known that, to capture the essence, 
one has to have a good physical intuition beforehand and 
the validity and accuracy of the wavefunctions strongly 
rely on this intuition and genius idea. From the first prin- 
ciples point of view, it is desired to get around this bias 
inherent in the variational approach as much as pos- 
sible to have a versatile method in hand. We need to 
construct wave functions which enhance the capability 
of removing biases posed on the variational forms. This 
is achieved at least partially by introducing enormous 
number of parameters. Recent development in the VMC 
method allows us to deal with a large number of param- 
eters.^^' Numerical techniques for optimizing 
many parameters are described in §4.2.5. By many pa- 
rameters in the part of refined single-particle wavefunc- 
tions as well as in the part of correlation factors such as 
Gutzwiller factor to punish two electrons on the same 
Wannier orbital, it opens a way of simulating strongly 
correlated systems by substantially reducing the biases 
as we see below. 

4.2.1 Functional form of variational wave functions 

The general functional form of wave functions we em- 
ploy is 

IV) = vm, (102) 

where \(j)) is a Hartree-Fock-Bogoliubov type wave func- 
tion called "one-body part," C is the quantum- number 
projector^^^' controlling symmetries of wave function, 
and V is the Gutzwiller- Jastrow factor^^^' includ- 
ing many-body correlations. In order to improve vari- 
ational wave functions within the sector classified by 
quantum numbers, we employ V that preserves symme- 
tries of jC\(j)). This means that C and V are commutable 
{VC = CV). 
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4-2.2 One-body part 

The one-body part usually corresponds to the mean- 
field Slater determinant with several variational param- 
eters. Though the Gutzwiller-Jastrow factor introduces 
many-body correlations through this variational hy- 
pothesis for the one-body part strongly restricts fiex- 
ibility of wave functions. The functional form of the one- 
body part has been reexamined and as many as possible 
variational parameters have been introduced in order to 
improve wave functions. 

Following this reexamination, we use a variational 
wave function in the form 



Impair) 



zeBZ 



feGAFBZ 



-| 

|o>, 

(103) 

where (p^^\k) and (p^'^\k) are k dependent variational 
parameters with the conditions 

^^^\-k) = ^^^\k) , ^^^\-k) = ^^^\k). (104) 

This form allows explicitly representing antiferromag- 
netic mean-field state with the periodicity Q by using 
the second term proportional to Lp^'^\k). Not only the 
Fermi sea state, this form also allows representing the 
BCS-type superconducting wavefunction with the pair- 
ing amplitude proportional to Lp^^\k). Therefore, dealing 
with Lp^^\k) and Lp^'^\k) directly as /c-dependent vari- 
ational parameters allows us to express various states 
such as paramagnetic metals, antiferromagnetically or- 
dered states, and superconducting states with any gap 
function within a single framework of Impair)- Moreover, 
since the number of the variational parameters increases 
scaled by the system size, it allows taking account of fiuc- 
t nation effects with short-ranged correlations. We call 
Impair) a "generalized pairing function" and Lp^^\k)., and 
Lp^'^\k) are called "pair orbitals." Introducing all the pos- 
sible ordered vectors Q would further generalize Impair)- 
However, this extension substantially increases the num- 
ber of variational parameters and computational costs 
(~ 0{N)). Therefore, one physically plausible Q has 
been attempted so far. 

By using the /c-dependent parameters Lp^^\k) and the 
Gutz wilier factor Vq defined below, our variational wave 
function can also represent the resonating valence bond 
(RVB) basis, which is known to offer highly accurate 
variational wave functions in spin systems. We note that 
the RVB basis can represent the state with spin correla- 
tions decaying with arbitrary power laws for increasing 
distance. 

For actual numerical calculations, we rewrite |(/)pair) in 
a real space representation: 



Impair) 



N/2 



|o) 



(105) 



with 



- y 



-iQ'Vj 



fcGAFBZ 



(106) 

Here, BZ means the summation in the Brillouin zone 
and AFBZ represents the folded zone for translational 
symmetry broken states with the periodicity vQ. One of 
the parameters {(^*^-^^(Aj), Lp^'^^k)} is not independent be- 
cause of the normalization of the wave function. We note 
that fij depends only on Vi — Vj because of the transla- 
tional symmetry. In practical calculations, the numbers 
of the variational parameters fij can be decreased to re- 
duce the computational load by restricting the range of 
fij into only short-ranged combinations. 

4.2.3 Gutzwiller-Jastrow factors 

In the variational study, the Gutzwiller-Jastrow type 
wave functions^^^' are often used to take account of 
many-body correlations. The Gutzwiller-Jastrow corre- 
lation factor V allows us to go beyond a single Slater 
determinant and a linear combination of many Slater de- 
terminants are generated after the operation of 7^ to a 
one-body wave function which is crucial in repre- 
senting strong correlation effects. Three types of many- 
body operators Vq^ ^d^i' ^J' which are called the 
Gutzwiller factor, the doublon-holon correlation factor, 
and the Jastrow factor, respectively, have been employed 
for the low-energy solvers so far. 

Gutzwiller has introduced a basic and efficient corre- 
lation factor Vg^^^^ which punishes double occupancy of 
up and down electrons at the same Wannier orbital as 



Vg = exp 



-9 



=n[i-a 



')nitni4 



(107) 

where ^ is a variational parameter and i represents the 
site and orbital indices. In the limit g ^ 00, Vg fully 
projects out the configurations with finite double occu- 
pancy as 



[1 - rii^riii 



(108) 



Vq is used for the Heisenberg model and the t-J model. 
In Hubbard- type models with finite U/t with the onsite 
interaction U and the typical transfer t, the double oc- 
cupancy at the same orbital and site is nonzero even in 
the insulating state. Thus we deal with Vg at a finite g. 

In order to take account of many-body effects beyond 
the Gutzwiller factor, the doublon-holon correlation fac- 
^Q^i33, 134 implemented in the wave function. This fac- 
tor comes from the idea that a doublon (doubly occupied 
site) and a holon (empty site) are bound in the insulator 
for large U/t. An example of the short-ranged correlation 
factor has a form 



Vd-h exp 



-iE4o\ 



(109) 
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where ai is a variational parameter. Here, ^^^^^ is written 

by 



(110) 



where the product runs over nearest-neighbor sites, 
and di = Ui^riii and hi = (1 — ni^){l — n^^) are dou- 
blon and holon operators, respectively. This factor takes 
into account the attraction of a doublon and a holon at 
the nearest neighbor site. The doublon- holon correlation 
factor Vd-h given by eq. (109) or by slightly different 
forms has been employed in several VMC studies. ^^^^^^ 
The correlation between doublons and holons at further 
distance may also be considered in more sophisticated 
forms.^^ 

Jastrow has introduced a long-ranged correlation fac- 
tor for continuum systems. This factor takes into ac- 
count correlation effects through two-body operators. In 
the Hubbard model at quarter filling, Yokoyama and 
Shiba have discussed the effects of the Jastrow-type cor- 
relation factor. Recently, Capello et al. have claimed a 
necessity of this factor to describe the Mott transition. 
The Jastrow factor Vj in lattice models has the following 
form: 



Vj = exp 



with two-body terms, where rii = rii^ + riii is a density 
operator and Vij = v{ri — rj) are variational parameters 
depending on the displacement — rj . 

In addition to the correlation factor, one can also oper- 
ate Hamiltonian matrix H to further approach the true 
ground state. This is the idea to implement Lanczos- 
type diagonalization partially. The first order correction 
is realized by operating 1 + aH with a variational pa- 
rameter a, which corresponds to the first order Lanczos 
step..80'i39 

4-2.4 Quantum-number projection 

In general, quantum many-body systems have several 
symmetries related to the Hamiltonian such as transla- 
tional symmetry, point group symmetry of lattice, U{1) 
gauge symmetry, and SU(2) spin-rotational symmetry. 
While symmetry breaking occurs in the thermodynamic 
limit, these symmetries must be preserved in finite many- 
body systems. 

Variational wave functions constructed from one-body 
parts and the Gutzwiller- Jastrow factors do not of- 
ten satisfy inherent symmetry properties, because the 
Hartree-Fock-Bogoliubov type one-body part comes from 
symmetry broken mean- field treatment. Even in the gen- 
eralized pairing wave function |(/)pair), the spin-rotational 
symmetry is broken by the orbital Lp^'^\k) which enables 
to include the mean-field AF state. 

The quantum-number projection technique^^^ enables 
to control symmetries of wave function. This technique 
has been used successfully in the PIRG method^^^ and 
the Gaussian-basis Monte Carlo method. ^^'^^ By us- 
ing the quantum- number projection together with the 
Gutzwiller- Jastrow factor, one can construct variational 



wave functions with controlled symmetries and many- 
body correlations. The quantum- number projection op- 
erator C is constructed by superposing transformation 
operators T^^^ with weights Wn'- 



£|(A) = 5]t^„T(")|0)=^ti;„|0(")), 



(112) 



where |(/)) and are the original one-body part and 

the transformed one-body parts, respectively. When C 
restores some continuous symmetry, the summation 
is replaced by the integration over some continuous vari- 
able. 

For example, the SU (2) spin-rotational symmetry pre- 
served in many effective models is restored by superpos- 
ing wave functions rotated in the spin space. The spin 
projection operator which filters out 5*^ = compo- 
nent of and generates a state with total spin S and 
= {) has a form 

j df2Ps{cosp)R{f2), (113) 

where i? = (a, /3, 7) is the Euler angle and the integration 
is performed over whole range of O. The weight Ps- (cos f3) 
is the ^-th Legendre polynomial. The rotational operator 
R{f2) is defined as 

(114) 



(111) R{n) = R'{a)Ry{l3)R'{j) = e'^^\'^^" e'^^^ 



where S'^ and are total spin operators of y and z 
directions, respectively. 

The one-body part introduced in this article contains 

only = component \^) = [E^J h4Ai\^^^\^y 
Then, the integration over 7 and a can be omitted and 
C^ld)) is written as 



= E 1^)^^ dl5smpPs{cosp){x\Rym<t>)- 

(115) 

with the complete basis set of the real space representa- 
tion \x). The integration over {3 is evaluated efficiently 
by the Gauss-Legendre quadrature in actual numerical 
calculations. Typically, for 6* = of the half- filled 
electron system for the single-band Hubbard model on 
square lattices with the sizes 4x4 and 14 x 14, 10 and 
20 mesh points are sufficient, respectively. 

There are also other quantum-number projections to 
restore symmetries. The total momentum projection 
and the lattice symmetry projection restore the transla- 
tional symmetry and the point group symmetry of lat- 
tice, respectively. The momentum projection is given by 
taking T^^^ as the operator to shift the state with the 
translation vector Rn. The state with the total momen- 
tum k is obtained by employing Wn = ex.p{ik • Rn) with 
the summation over n in eq.(113) for all the possible spa- 
tial translations in the finite size system. The momentum 
projections can be redundant with the spin projection, if 
the translational symmetry is automatically restored by 
the spin projection, where a superposition of the spin- 
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rotated wave functions includes a superposition of spa- 
tially translated wave functions. In such cases, the mo- 
mentum projection does not make a better wavefunction 
any more and can be omitted. 

With the quantum number projections, not only the 
ground state but also excited states with specified quan- 
tum numbers are obtained. It is useful, for instance, in 
obtaining the dispersion of quasiparticles by taking the 
momentum projection. 

4-2.5 Optimization method 

Here, we summarize the basic idea of wave func- 
tion optimizations. The stochastic reconfiguration (SR) 
method^^^ introduced by Sorella has been employed in 
many-parameter optimization. 

[1] Wave function optimization by energy mini- 
mization 

We first recollect the conventional way of minimizing the 
variational energy = {i^a\y^\i^cx) / (i^cxli^a.) estimated 
from the wave function \tjjot) with variational parameters 
{ak\k = 1, • • • Here a denotes the initial vector in 
the p-dimensional parameter space. 

The energy Eoc+^ is expanded up to the second order 
around a: 



-E 



9klk- 



1 ^ 

is 



hM-ikli^O{-i^), (116) 



where 7 is the vector for parameter variations, 
d 



9k 



dak 



Eoc (/c = 1, • • • ,p) 



are the energy gradient vector gf, and 



h 



ki 



dai 



E^ {k,£ = ,p) 



(117) 



(118) 



are the elements of the energy Hessian matrix h. 

The steepest decent (SD) method gives the updated 
variational parameter by 

<^fc = Q^fe + 7fc, (119) 
where the change from the initial value ak should be 

% = -i2^M9i {-f = -X-'g). (120) 

with a suitably chosen nonsingular matrix X. This gen- 
eral form reduces to the steepest descent (SD) method 
in the choice Xk£ = Ski/^t and to the Newton method 
by taking the hessian for X as Xk£ = hk£^ respectively. 
In general, X should be properly chosen to accelerate the 
optimization. 

[2] Stochastic reconfiguration method 

Sorella has introduced the SR method^^^ for a better 
choice of X in optimizing many variational parameters. 
For the normalized wave function 

\^c.)= (121) 



up to the first order around a, IV^^+'y) is expanded as 

p 



where |V^/cq-) (/c = 1, • • • ,p) are the derivatives of 

1 f_d_ 



koc) 



(123) 



The wave function set {\tlJka)\k = I,-- - ,p} forms 
nonorthogonal basis in the p-dimensional parameter 
space. The norm change from |?/;q,) to IV^a+'y) is 

2 



^ Ikliii^kcxli^icx) = ^ IkliSki' (124) 



kl^l 



kl^l 



Equation (124) shows that S whose {k^i) element is Sm 
is the metric matrix in the parameter space. 

The SR method chooses S as the matrix X in eq. (120), 
namely 



7/c 



-At 



p 



(125) 



where At is a small constant. The SR method takes into 
account the variation of the wave function in addition 
to the SD method. We can derive eq. (125) by mini- 
mizing the functional J^sR = ^^^lin. + ^^norm with a 
Lagrange multiplier A. Here Z\£^iin. = ^^Qklk is the 
linear change of the energy. The stationary condition 
dJ~SR/djk = (^ = • ' ' 7^) leads to the SR formula 
(125) with At = (2A)"^ The SR method is more sta- 
ble than the conventional method, because the SD and 
hessian methods sometimes cause a large change in the 
wavefunction even though the changes in the variational 
parameters are small. This large change in the wavefunc- 
tion causes an instability in the iteration, which requires 
to keep At very small and the iteration becomes inef- 
ficient. The SR method solves this difficulty. To avoid 
the numerical instability possibly caused by an extremely 
large in eq.(125), it is also useful to take (1 + s)Skk 
instead of Skk for its expression in eq.(125) with a small 
constant e}^^' 

The positive definite matrix of S may have an eigenvec- 
tor with very small eigenvalues after the diagonalization. 
The variation in the direction in such an eigenvector is 
redundant and may be truncated for a better efficiency. 
Typically, the parameters are taken as At = 0.1, e = 0.2, 
and ^wf = 0.001. For more details including the practical 
implementation of the optimization, readers are referred 
to Ref.80. 

4.2.6 Benchmark 

The accuracy of the multi- variable VMC method has 
been critically tested in various cases and in many cases 
has proven to be a good low-energy solver as comparably 
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Fig. 16. (Color online) Total energy E/Ns as a function of 1/L^ 
for the Hubbard model on a square lattice with Ns = L x L lat- 
tice with the transfer t = 1 to the nearest neighbor only, and the 
onsite interaction U = 4 at half filling n = 1.^^ The exact values 
are calculated by the exact diagonalization (L = 4) and AFMC 
(L = 6,8,10,12,14). Error bars are comparable to the symbol 
size. The accuracy is enhanced with improving variational forms. 
The best results by (purple) square are obtained by operating 
the spin quantum number projection to = 0, Gutzwiller pro- 
jection, doublon-holon correlation factor, Jastrow factor and the 
1st order Lanczos step to the paired singlet function \(f)pair)- In 
this case, the relative error is smaller than 0.5 % irrespective of 
the system size. 




5=0 S=l S=2 S=3 

Total Spin 

Fig. 17. (Color online) Excitation spectra of ID Hubbard model 
compared with the exact diagonalization result (ED). VMC can 
reproduce the lowest energy states of each quantum number. 
VMC+LS represents the results obtained by the form corre- 
sponding to the squares in Fig. 16. Energies are sufficiently accu- 
rate without the Lanczos step as we see for circles. 



accurate as PIRG reviewed below. In Fig. 16, we show the 
energy accuracy for the case of the Hubbard model on 
the square lattice with various system sizes at the onsite 
interaction U = 4 and the nearest-neighbor transfer t = 
1.^^ Excitations are also tested in Fig. 17, where different 
total spin momentum and parity states are obtained 
by the present VMC for the Hubbard model on a ID 8 x 1 
lattice at t/ = 4 and t = 1. 

4.3 Path-integral renormalization group 

Now we introduce PIRG method applied for low- 
energy effective models. Since reviews are also found in 
the literature, we review only the essential part here. 
This method allows to approach the ground state of the 
model with a high accuracy. It starts, as in other Monte 
Carlo and projection methods, with a relation 

\^g) = lim e-"^|(/)initiai). (126) 

r^oo 

for an arbitrary chosen state |(/)initiai) that is not 
orthogonal to the ground state li^g). By follow- 
ing the Feynman path integral, the operation of 
exp[— ri^] is decomposed with Ar as exp[— ri^] ~ 
[exp[— Ari^x] Yli exp [-ArHu,]]^. Here we take a suf- 
ficiently large N so as to satisfy r = AfAr. For sim- 
plicity, we have taken an example of the Hubbard type 
model with the onsite interaction at the ith site. 
When we take a Slater determinant | ^initial), the oper- 
ation of exp[— Ari^x] to the Slater determinant simply 
generates another single Slater determinant. However, if 
we operate exp[— Ari^f/.] to a single Slater determinant, 
the result is a linear combination of two Slater determi- 



nants when this interaction is transformed by the discrete 
Stratonovich-Hubbard transformation. To approach 
the ground state we need to operate exp[— Ari^f/.] many 
times, which requires a linear combination of an expo- 
nentially large number of Slater determinants. It easily 
exceeds the accessibility by computers. Then by restrict- 
ing within the computationally tractable range, a linear 
combination of L Slater determinants in a partial Hilbert 
space truncated from the original Hilbert space as 

|^(^)) = X:c«|</.i^)>, (127) 

is employed in PIRG. The coefficients Ca and the choice 
of nonorthogonal basis ^i^'^ are numerically optimized. 
Then by increasing L systematically, the ground state is 
speculated from the extrapolation to large L (essentially 
the hmit / 00) When the Hilbert space is truncated, it 
can lose the original symmetries of the Hamiltonian that 
are guaranteed by the conservation law, while the ground 
state should be an eigenstate of good quantum numbers 
that are conserved in the original Hamiltonian. For in- 
stance the total spin and total momentum are normally 
good quantum numbers in the Hubbard-type Hamilto- 
nian and in the ground state, one of the quantum num- 
bers has to be chosen for each conservation law. To re- 
store the original symmetry, quantum number projection 
may be imposed to keep the quantum number of the 
truncated state. This quantum number projection was 
combined with PIRG that has proven much better accu- 
racy. Since PIRG does not suffer from the negative sign 
problem known in AFMC, a high- accuracy calculations 
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have been made possible. We will refer to applications to 
ab initio calculations in §5. 

5. Applications 

Applications of the three-stage RMS formalisms com- 
bining ab initio electronic structure calculations, down- 
folding and low-energy solvers are diverse and it is not 
possible to cover all of them in this review. Here, we just 
pick up several examples to demonstrate the efficiency, 
accuracy and versatility of the method in practical ap- 
plications. 

5. 1 Dynamics of semiconductors 

Despite its success in sp semiconductors and insula- 
tors, the LDA is not satisfactory when it comes to elec- 
tron excited states. The LDA band gap is 1/2^2/3 of 
measured values in a wide range of materials. The cause 
of the deviation is not ascribed to approximations in the 
exchange correlation functional. In fact, the GGA yields 
similar results as the LDA and does not cure the band 
gap problem. The problem originates from the fact that 
Kohn-Sham eigenvalues do not correspond to observable 
quantities. Schematically it is expressed as Fig. 18. The 
Kohn-Sham eigenvalues represent energy levels of a 
electron system, where the lowest N levels are occupied. 
On the other hand, what is observed in (inverse) photoe- 
mission measurements is electron addition/removal en- 
ergy, which is the total energy difference between the N 
electron system and the ± 1 electron system. Theoret- 
ically, they are obtained by the spectral function of the 
one particle Green's function, which can be computed 
from first-principles, e.g. in the GW approximation. 
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Table I. Band gap of silicon and diamond in LDA, Hartree-Fock, 
Hartree and GW approximations. The negative values mean the 
band overlap in metals. The energies are given in units of eV. 

term is replaced with the screened one, thereby the in- 
crease of the gap is suppressed. As a result, the GW 
yields the band gap close to the experimental one (Table 
I). The same trend is observed also in diamond. 

The GW self-energy is nonlocal and energy dependent. 
The nonlocality increases the band gap substantially and 
solves the underestimation of the gap in LDA. The en- 
ergy dependence, on the other hand, reduces the gap and 
partially cancels the nonlocal effect. Spatial and energy 
dependence of the dielectric function e(r, r'; u) is also im- 
portant. If the dielectric function is approximated as a 
function of |r — r'|, the value of the gap changes signifi- 
cantly (local field effect). The energy dependence cannot 
be neglected either, but a simplification using the plas- 
mon pole model is valid in weakly correlated materials. 
In the plasmon model, the imaginary part of the dielec- 
tric function is approximated as a single delta function. 
This enables us to extrapolate the energy dependence of 
the dielectric function from its static value. Many GW 
approximations adopt this approximation in order to im- 
prove computational efficiency. 



(a) 




O 



(c) 



Fig. 18. (a) Kohn-Sham energy is the energy level of a non- 
interacting electron system, where the lowest N states are occu- 
pied, (b) Quasiparticle energy is the electron removal / addition 
energy to the ground state of the N electron system, which is 
different from (a) when electron interaction is considered. Opti- 
cally excited state is shown in (c). Both an electron and a hole 
exist, and their interaction brings another many-body effect. 



The band gap of silicon is, for example, 1.17 eV exper- 
imentally. The LDA gap is 0.5 eV, much smaller than the 
experiment. In the Hartree-Fock approximation, the gap 
increases substantially to > 6 eV because of the nonlocal 
exchange term. In the GW approximation, the exchange 
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Fig. 19. (Color online) Band gap of semiconductors and insula- 
tors. Tj^g dotted line is the ideal line that the theoretical 
gap agrees with the experiment. The LDA gap is systematically 
smaller than the experimental one, which is greatly improved by 
adding the GW self-energy correction. 



After a seminal work by Hybertsen and Louie, many 
groups applied the GW method to a number of materials. 
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including solids, surfaces, molecules, clusters and nanos- 
tructures. By today the GW method is established as a 
reliable method to predict the band gap of weakly and 
moderately correlated materials within 10-15 % error to 
experimental values (Fig. 19^^^). 

On the other hand, the conventional GW method is 
not satisfactory for strongly correlated materials. This is 
partly because the DFT gives a poor starting Hamilto- 
nian. The band gaps are too small, and sometimes in- 
sulators are wrongly described to be metallic. Moreover, 
localized electron levels are too shallow in many cases. 
These drawbacks result in over-screened Coulomb inter- 
action W, which in turn yields inaccurate self-energy. 
In principle, the dressed Green's function (including the 
GW self-energy correction) can be used to recalculate 
the self-energy. The calculation is continued until the 
self-consistency is achieved. The fully self-consistent GW 
calculations have not been performed extensively up to 
now, partly because they are both computationally and 
technically demanding. The fully self-consistent GW cal- 
culations for the electron gas, however, have been per- 
formed in detail with a rather discouraging result with 
regard to the excitation spectrum. On the other hand, 
the total energy is found to be in almost perfect agree- 
ment with the quantum Monte Carlo. Applications to 
silicon have been performed with four different methods, 
but the results are inconsistent. ^^^^^^ Further system- 
atic calculations are anticipated for deeper insight into 
self-consistency. 

Another approach to reach self-consistency is to use 
the GW self-energy to update the one-particle wavefunc- 
tions and eigenvalues. Due to the energy dependence of 
the self-energy it is not clear, however, how this should be 
done. Consider the self-energy in the Kohn-Sham basis, 

(V^kn|S(a;)|V^km) . (128) 

To use this self-energy in place of the LDA exchange- 
correlation potential, it is necessary to determine the en- 
ergy in some way. We note that the self-energy matrix 
is required to be Hermitian in order to produce a set 
of ort honor mal wavefunctions. There are several choices. 
The simplest one is to fix the energy at some chosen en- 
ergy, say, the Fermi energy or the center of the band of 
interest. Another choice is to take the average^^^ 

(V^knlKclV^km) = ^ [(V^kn|S(ekn)|V^km) 

+ (V^knlS(ekm)lV^km)]. (129) 

This scheme is called quasiparticle self-consistent GW 
(QSGW) method. The scheme was applied to many 
materials including transition metal mono-oxides and / 
electron systems, and improvement over LDA was con- 
firmed in strongly correlated materials. A self- 
consistent scheme based on the static COHSEX approxi- 
mation (eq.(34)) is also reported. Recently, Sakuma et 
al proposed an alternative scheme. In their scheme, 
the quasi-particle equation is solved with neglecting the 
imaginary part of the self-energy, 

det \uj - Ho{k) - J?E(k,cj)| = . (130) 
The obtained quasiparticle wavefunctions {i^^^} 
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orthonormal because of energy dependence of the self- 
energy. The wavefunctions are then orthonormalized us- 
ing Lowding's scheme, that generates the closest set 
of orthonormal orbit als {V^kn} to {V^kJ}' where 

= jP^^C , (131) 

CC^ = 5 . (132) 

Smn (^kml^kn) • (133) 




(^'^) Energy [eV] 

Fig. 20. (Color online) Optical absorption spectra of (a) silicon 
and (b) LiF.^ The dashed line is the calculated spectrum in the 
independent particle approximation using Kohn-Sham wavefunc- 
tions and GW quasiparticle energies, the solid line includes the 
electron-hole interaction in the GW-BSE scheme, and open and 
closed circles are experimental data. The BSE spectrum is red 
shifted compared to the independent particle approximation due 
to electron-hole interaction, and agrees well with the measure- 
ment. 



Optical absorption spectrum is another quantity where 
electron interaction is crucial both quantitatively and 
qualitatively. Figure 20(a) shows the optical absorption 
spectrum of crystalline silicon. The dashed line is the 
spectrum in the independent particle approximation, 

/ D \ 2 occ unocc 
^ ^ V c 

xS{uj-{Ec-E,)), (134) 

where Ey and Ec are the GW quasiparticle energies. As 
the GW gap is accurate, the threshold of the spectrum is 
close to the experimental one. However, the peak position 
is too high compared to the experiment. Moreover, the 
first peak at 3.5 eV is not reproduced in the calculation. 
These discrepancies clearly show that many-body effects 
are crucial for the optical absorption spectrum. 

The spectrum including many-body effects is obtained 
by the imaginary part of the macroscopic dielectric func- 
tion, 

eM(^) = lim — r • (135) 

The independent particle approximation eq.(134) adopts 
two approximations on top of eq.(135). One is that the 
G 7^ components are neglected when e is inverted. 
Namely, the local field effect is neglected in the calcula- 
tion. Analysis revealed that this effect is minor in sili- 
con.^^^ The other approximation is the RPA for the po- 
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larization eq.(31). To go beyond the RPA and include 
electron-hole interaction, we start with Hedin's equa- 
tion. Firstly, the vertex function is evaluated in the fol- 
lowing way. The GW approximation is adopted in the 
second term of eq.(25). Then it follows that ST^/SG = 
iW + iG{SW/SG). Assuming that the screening effect is 
not affected by electron- hole excitations, the SW/SG con- 
tribution is safely neglected. Putting thus obtained ver- 
tex function into eq.(24), we can compute an improved 
polarizability. The final form called Bethe-Salpeter equa- 
tion is written as""^^""^' ""^^^ 



lim [vG^o{q)PG=G' 



.o(q,c^)], (136) 



P(1,2) = 4P(1,1,2,2), 



(137) 



^P(l,2,3,4) = ^Po(l, 2,3,4) ■ 



^Po(l,2,5,6) 



X K{5, 6, 7, 8)'^P(7, 8, 3, 4)d(5678), (138) 

which is diagrammatically illustrated in Fig. 21. Here 
"^P is a correlation function between the electron and the 
hole. The first term in eq.(138) is the correlation function 
in the independent particle approximation, whereas the 
second term represents the electron-hole interaction. It is 
characterized by the electron-hole interaction kernel 
given by 



6{l,2)5{3,4-)v{l,3) 
5(1, 3)5(2, 4) 1^(1+, 2), 



(139) 



where ^ is a modified bare Coulomb interaction, in which 
G = component in the Fourier representation is re- 
placed with 0. The second term, called direct term, rep- 
resents electron- hole attraction (excitonic effect), while 
the first term (exchange term) is the local field effect. 
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Fig. 21. (Color online) Feynman diagrams representing Bethe- 
Salpeter equation. The second term in the right hand side is the 
electron-hole interaction (top panel). The interaction is charac- 
terized by the kernel (bottom panel), where the second term 
represents the excitonic effect. 



Further approximations are introduced to make first- 
principles calculations to real materials feasible. The ex- 
cited state is approximated to be a linear combination of 
single electron-hole excitations (Tamm-Dancoff approxi- 
mation), 

elec hole 

\S)=J2J2^cAbl\0), (140) 



where aj and bl are creation operators for an electron and 
a hole, respectively. In addition, the frequency depen- 
dence of the screened Coulomb interaction is neglected 
and its static value is used. With these approximations, 
the Bethe-Salpeter equation (138) is reduced to the fol- 
lowing eigenvalue problem. 



{Er-Er)At 



S 

c'v' 



(141) 



^/'e*(ri)^/'.(r2)if(l,2,3,4) 



X t/-e'(r3)t/':,(r4)d(1234). 



(142) 



Using the eigenstates and eigenvalues of eq.(141), the 
macroscopic dielectric function is obtained by 

62(c.) = f — E KOI^A • ^m^Kuj - Qs) . (143) 
\mujj ^ 

Many-body theory for optical absorption is seen al- 
ready in 1960's.^^^ In the 70's, semi-quantitative calcu- 
lation was reported by Hanke.^^^ First ab-initio calcula- 
tion was carried out in 1995 for sodium cluster. The 
scheme was applied to solids 2-3 years after. ''^ 

The Bethe-Salpeter-Equation (BSE) spectrum of Si is 
shown by the solid line in Fig. 20 (a). The first peak ap- 
pears and the second peak is red shifted by the electron- 
hole interaction. The spectrum compares quite well with 
the experiment. Good agreement with experiment is also 
seen in an insulator with a large band gap. Figure 20(b) 
shows the spectra for LiF. In sharp contrast to silicon, a 
bound exciton peak is formed in the gap by the inclusion 
of electron-hole interaction. Consequently, the threshold 
of the spectrum is red shifted substantially.^ 

A few different approaches have also been developed 
for the optical absorption spectra. Time Dependent Den- 
sity Functional Theory (TDDFT)^^^ is studied inten- 
sively in the last decade. It was pointed out by Runge 
and Gross that by taking the time- dependent one elec- 
tron density as a basic variable, DFT can be rigor- 
ously extended to treat dynamic response of many elec- 
tron systems. The many-body effects are included in the 
exchange-correlation kernel, which is the density deriva- 
tive of the exchange correlation potential. Most applica- 
tions to real materials adopt the adiabatic local density 
approximation (ALDA) that neglects the nonlocal effects 
in both time and space. The TDDFT using the ALDA 
works well for finite systems. As the system size becomes 
larger, discrepancy from the experiment becomes signif- 
icant. Improvement of the exchange-correlation kernel is 
needed for application to solids. 

Nakamura et al. followed the three-stage scheme for 
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the optical absorption spectrum of GaAs.^^'' Start- 
ing with the GGA band structure, they derived a 
low-energy tight-binding Hamiltonian using the maxi- 
mally localized Wannier function procedure. The on- 
site Coulomb interaction is evaluated by the constrained 
DFT method. Solving the derived model using the 
Hartree-Fock approximation supplemented by the single- 
excitation configuration-interaction method considering 
electron-hole interactions, they obtained the spectrum in 
good agreement with experiments. 

5.2 3d transition metal and its oxides 

5.2.1 Transition metal 

Simple substances of the 3d transition metals are rea- 
sonably described by the LDA. They are classified as 
moderately correlated materials. The value of U is esti- 
mated to be 3-5 eV in the cRPA,^^' which is consistent 
with the above picture. 

One of the problems in the LDA band structure is too 
wide d band width. A GW calculation of Ni by Aryaseti- 
awan^^ showed that the GW self-energy correction raise 
the bottom of the d band by about 1 eV, resulting in a 
band narrowing, in agreement with experiments. How- 
ever, experimentally observed satellite at -6 eV below 
the Fermi level is not reproduced even in the GW. In 
addition, the exchange splitting does not change signifi- 
cantly compared to LDA, and larger than measured val- 
ues. In fact, the satellite originates from short-range cor- 
relations, which is not properly described in the GW ap- 
proximation. The problem was solved by a T-matrix cal- 
culation.-*^^ The problem of the exchange splitting, as well 
as the 6 eV satellite and band narrowing, was settled 
down later by the LDA+DMFT,i^^ and GW+DMFT 
calculation. 

Another well-known problem in the LDA is mag- 
netism. In iron, for example, the nonmagnetic hep struc- 
ture becomes more stable than the ferromagnetic bcc 
in the LDA. The GGA correctly describes the ground 
state in this particular case,^^^^*"^ but in general care- 
ful analysis is needed for magnetic properties. Looking 
at finite-temperature properties, one needs a formalism 
that takes into account the existence of local magnetic 
moments above the Curie temperature Tc. Many-body 
effects which incorporate the local atomic character of 
the electrons are essential. This can be achieved by the 
LDA+DMFT method. Applications to Fe and Ni repro- 
duced semi-quantitatively the ferromagnetic susceptibil- 
ity above Tc and temperature dependence of the ordered 
moment below Tc.^^^ 

5.2.2 SrVOs 

Transition metal perovskite compounds exhibit var- 
ious intriguing electronic and magnetic properties. 
Among them, SrVOs can be regarded as a prototype. 
The material is cubic. There is no GdFeOa-type lattice 
distortion, and only one formula unit is contained in the 
unit cell. 

Experimentally, the compound is a paramagnetic 
metal. Fujimori et al. found by photoemission spec- 
troscopy (PES) that the occupied d band has a double- 
peak structure, one within about 1 eV of the Fermi level 




R r x M r 

Fig. 22. LDA band structure of SrVOs. The three states crossing 
the Fermi level have strong V-t2p component. There are nine 
states at [-8eV:-2eV], which is mainly of oxygen 2p character. 

with a sharp Fermi cutoff, whereas the other centered 
at ~1.5 eV below the Fermi level. ^''^ This suggests that 
SrVOs is a correlated metal, where the peak around the 
Fermi level corresponds to the quasiparticle peak, and 
the other peak is the lower Hubbard band. The inverse 
photoemission spectrum shows a peak at 2.5-3 eV, which 
can be interpreted as an upper Hubbard band.^^^ 

Fujimori et al. also studied other d^ electron sys- 
tems including ReOa, VO2, SrVOs, LaTiOa, YTiOa, and 
found that as the bandwidth decreases, deviation from 
the band structure calculation becomes substantial, and 
the weight near the Fermi level is transferred to the 
higher binding energy. 

Comparison between SrVOs and CaVOa has also at- 
tracted much attention. Both compounds are metallic 
perovskites having a single d electron. A noticeable dif- 
ference is that the CaVOa has a distorted structure. This 
leads to reduction in hopping between the t2g orbitals of 
neighboring V atoms mediated by the 0-p orbital, which 
would enhance the correlation effects. Early PES experi- 
ment for Cai-a^Sra^VOa reported that the spectral weight 
is indeed transferred from a coherent to incoherent peak 
as X decreases. However, later high-energy photoemis- 
sion experiment casted doubt on this conclusion, finding 
that the bulk spectra is insensitive to x.^^^ A more recent 
low-energy PES measurement supported this result. 
In the measurement, it was also found that the spec- 
tral intensity is suppressed near the Fermi level. This is 
consistent with the cluster extension of DMFT for the 
Hubbard model showing the suppression of the coherent 
peak with a pseudogap formation as we discussed and 
illustrated in Fig. 15).^°^ 

In the LDA band structure (Fig. 22), three states hav- 
ing strong y-t2g character cross the Fermi level. They 
form an isolated band, below which is an oxygen band 
at [-8eV:-2eV], while the V-e^ band is just above the 
t2g band. Two valence electrons are transferred from V 
to oxygen, thus a single electron occupies the t2g band. 
The LDA correctly reproduces paramagnetic metal na- 
ture. However, the t2g band is much broader than the 
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experimental quasiparticle peak: The latter is about 60 
% of the former in width. In addition, the satellite struc- 
ture at -1.5 eV does not exist in LDA. LDA+U cannot 
reproduce these features neither. A standard 1-shot GW 
reduces the band width by about 30 %, in reasonable 
agreement with the experiment. 

Coexistence of coherent and incoherent peaks can be 
reproduced only by going beyond static mean-field treat- 
ment of electron correlation effects. DMFT is a possible 
solution for this. Liebsch performed a DMFT calculation 
of SrVOs and CaVOs using the tight-binding Hamilto- 
nian fitted to the LDA t2g bands. For a reasonable 
choice of U ^ 4eV, a narrowing of quasiparticle peak 
and evolution of Hubbard bands are observed, in agree- 
ment with experiments. It was also found that the or- 
thorhombic distortion causes a weak transfer of spectral 
weight from the coherent to the incoherent peak. Later 
on, Pavarini et at. studied SrVOs and other three per- 
ovskites, CaVOa, LaTiOs and YTiOs, ranging from cor- 
related metal to magnetic insulator. They first extracted 
t2g bands using NMTO-Wannier procedure.^ ''^ Then the 
derived low energy Hamiltonian, with several values of 
U between 3-6 eV, was solved by DMFT. They found 
that the main features of the photoemission spectra for 
all four materials, as well as the correct values of the 
Mott-Hubbard gap for the insulators were reproduced by 
taking U to be 5 eV. Both SrVOa and CaVOs are corre- 
lated metals, while the quasiparticle peak disappears and 
the system becomes insulating in LaTiOs and YTiOa, as 
can be seen in Fig. 23. It was also revealed that the lat- 
tice distortion leads to reduction of not only band width 
but also of effective orbital degeneracy, which plays an 
important role in the metal-insulator transition. 

For a full quantitative treatment, ab-initio determi- 
nation of U is important. The value is estimated to be 
3.0-3.5 eV in cRPA,^^'^^ which is smaller than that used 
in previous LDA+DMFT calculations. This discrepancy 
implies that spatial correlations significantly enhance the 
electron correlation effects and tendency for the Mott 
insulator. Careful ah initio analysis on e.g., long-range 
interaction and non-local self-energy effects, is an open 
issue. 

5.2.3 VO2 

Vanadium dioxide is a material under debate for many 
years. As the temperature decreases, the material shows 
metal-insulator transition at 340 K from the high tem- 
perature metallic phase in the rutile structure to the low 
temperature insulating phase with the monoclinic (Ml) 
structure. There has been long discussion about the 
transition, with particular interest on the role of electron 
correlations in forming a gap. 

Seen from the band picture, the electronic states are 
understood as Fig. 24.^^^ The low-energy states near the 
Fermi level are of strong vanadium 3d character. The 
crystal field makes the d states split into t2g and Cg. Since 
the structure is not cubic, the t2g states are lifted further 
into Cg and aig state. The isolated vanadium atom has 
three d electrons. Two of them are transferred to oxy- 
gen 2p orbitals in VO2, thereby VO2 is a d^ system. The 
remaining d electron partially occupies the aig band so 
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Fig. 23. (Color online) (a) Occupied Wannier orbital of LaTiOs 
in primitive cells (right) and a subcell (left) obtained by 
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Fig. 24. Orbital energy diagram VO2 in monoclinic phase. The 
degeneracy of the V-t2^ states are lifted by the non-cubic crystal 
symmetry. The Peierls distortion leads to a coupling between two 
a\g states each on the neighboring V sites, forming the bonding 
and the anti-bonding states. The bonding state accommodates 
two electrons, and the system becomes gapful. 



that the system is metallic in the rutile phase. In the Ml 
phase, two vanadium atoms form a dimer. This Peierls 
distortion causes strong hybridization between the a\g 
orbitals of the two vanadium atoms. Then the bond- 
ing state is fully filled, which opens a gap between the 
bonding a\g and unoccupied band. The overall fea- 
ture was confirmed by first-principles calculations^^^' 
in the LDA. However, it is also found that the bond- 
ing a\g band overlaps with the band, yielding metal- 
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lie behavior in eontrast with the experiment. This may 
be ascribed to the band gap problem of LDA. If we in- 
clude many-body effects, the may shift up and the 
gap would open. On the other hand, there is another 
phase in which one half of vanadium atoms dimerize, 
while the other half form chains with equal space. This 
M2 phase is also insulating, which suggests that VO2 
may be a Mott insulator. Some authors claimed that the 
electron correlation plays a major role.^^"^' The con- 
troversy is not yet settled down and correlation effects 
beyond LDA are discussed with various techniques such 
as LDA+DMFT,i^6 simplified GW scheme^^^ or full 
ab-initio GW.i^O'i^i 

Here we show how the GW works. Figure 25(a) shows 
electronic structure of Ml phase obtained by LDA. As 
described above, the bonding aig orbitals are near the 
Fermi level, which are located in [-0.5 eV: eV]. These 
states are almost fully filled, but there is an overlap with 
the Cg band that is located just above the aig. As a result, 
there is a small hole (electron) pocket in aig (ep band, 
and the system becomes metallic in LDA. 

Now we add the self-energy correction to the LDA 
Kohn-Sham energies. Figure 26 shows the self-energy as 
a function of energy for selected states. The self-energy 
does not decrease monotonically but has dips and peaks. 
This behavior is quite different from weakly correlated 
semiconductors, such as silicon. Because of the pecu- 
liar energy dependence we need to treat full energy de- 
pendence of the self-energy. In fact if we compute the 
GW quasiparticle band by linearizing the energy depen- 
dence of the self-energy, as most ab-initio GW calcula- 
tions assume, the aig band gets too narrow. Also the 
non-linearity yields a weak satellite structure in the one 
electron spectral function above the Fermi level, but not 
below, in contrast with the LDA+DMFT result. 

Quasiparticle band structure is plotted in Fig. 25(c) by 
circles. We can see that band overlap between the aig and 
Cg is removed by the self-energy correction. It should be 
noted that the Kohn-Sham wavefunctions have too much 
hybridization between aig and near the A point, there- 
fore off-diagonal elements are essential to disentangle the 
bands. If we neglect the off-diagonal self-energy elements, 
the aig and overlap each other even in the GW level, 
as shown in Fig. 25(b). The quasiparticle aig band is now 
isolated and direct gap opens. For the opening of the 
fundamental gap, the effect of self-consistency plays a 
crucial role. 

Biermann et al. studied the compound by LDA com- 
bined with a cluster extension of DMFT.^^^ Starting 
from the LDA band structure, they extracted three t2g 
states per V atom, and constructed a multi-band Hub- 
bard model. The model was solved using cluster DMFT 
including all off-diagonal terms in orbital space. More 
precisely, instead of calculating the self-energy from a lo- 
cal impurity model embedding one single atom in a self- 
consistent bath, a pair of V atoms in a bath is considered. 
This is important because the formation of singlet pairs 
resulting from the strong dimerization can be captured 
only in a cluster extension. 

They carried out calculations for both rutile phase 
and Ml phase. The calculated spectral function is shown 




Fig. 25. (Color online) Electronic structure of VO2 in the mono- 
clinic phase, (a) The LDA band shows metallic character, (b) The 
GW self-energy correction with diagonal elements (in the Kohn- 
Sham basis) pulls up the conduction band, while the bands are 
still entangled at around the A point, (c) The bands are disen- 
tangled by including the off-diagonal elements, consequently the 
bonding aig band is isolated. ^^-"^ 
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Fig. 26. (Color online) GW self-energy of VO2 in the monoclinic 
phase at F point in the aig band.-*^^-*^ The self-energy is not a 
smooth function of frequency, which yields a satellite structure 
in the spectral function. The straight line with a positive slope 
is a; — e^^"^. The intersection between the line and the real part 
of the self-energy gives the quasiparticle energy. 
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Fig. 27. Spectral function of t2g orbit als for rutile and Ml phases 
obtained by the LDA+DMFT calculation^^^ The lower Hubbard 
band is located at -1.2 eV in the R phase, while the prominent 
peak at -0.8 eV in the Ml phase is the quasiparticle peak. 

in Fig. 27. For the rutile phase, the results of single-site 
and cluster-DMFT calculations are very similar. A clear 
quasiparticle peak is found near the Fermi level, and 
many-body effects reduce the bandwidth. Hence, the ru- 
tile phase can be characterized as a metal with inter- 
mediate correlation. Hubbard satellites are observed at 
high energies at -1.5 eV below and 2.5-3 eV above the 
Fermi level. In the Ml phase, nonlocal self-energy opens 
up a gap of about 0.6 eV (for U =4 eV and J=0.68 eV), in 
reasonable agreement with experiments. There is a sharp 
coherent peak at -0.8 eV. Below this peak is a weak lower 
Hubbard satellite at -1.8 eV, whereas the broad peak cen- 
tered at 2.2 eV is the upper Hubbard band. Charge dis- 
tribution is modified significantly, and the single electron 
occupies almost entirely the aig orbital. The low-energy 
nature of the insulator is quite different from that of 
a standard Mott insulator in which local moments are 
formed. In fact, at low frequency, the onsite component 
of the self-energy for the aig orbital behaves linearly as a 
function of frequency in contrast to the l/u behavior for 
the local moment Mott insulator. Based on these results, 
they concluded that at low energy, the compound is a 
Peierls insulator assisted by strong Coulomb correlation. 

The gap in the M2 phase appears to be a correlation- 
origin Mott gap, but it is not a settled issue. 

5.24 Sr2V04 

Sr2V04 has a layered perovskite structure and is iso- 
morphic to the mother compounds of a cuprate supercon- 
ductor La2Cu04.^^^' This compound has one 3d elec- 
tron per V site {d^ system) with strong two-dimensional 
anisotropy and has a dual relation to the one 3d hole 
per Cu sites {d^ system) of the cuprates. The duality 
is, however, not perfect because, in Sr2V04, the orbital 
degeneracy of d^ electron remains between dyz and dzx 
orbitals. The crystal field splitting of dxy orbital is also 
rather small (~ 0.08 eV in the LDA calculation), which 
evokes us importance of orbital physics. 

Since the 3d t2g bands are located near the Fermi level 
and are rather isolated from other bands as we see in 
Fig. 28, a low-energy effective model of the form (71) 
for the t2g Wannier orbitals has been derived. Af- 
ter the downfolding, the onsite interactions among the 



Wannier orbitals of intraorbital xy^yz{zx) and interor- 
bital xy-yz{xy-zx) and yz-zx combinations are U = 
2.77,2.58,1.35 and 1.28 eV, respectively. The onsite ex- 
change interactions between xy-yz{xy-zx) and yz-zx or- 
bitals are 0.65 and 0.64 eV, respectively. The nearest 
neighbor transfers between xy-xy, yz-yz and zx-zx or- 
bitals in X direction are -0.22, -0.05 and -0.19 eV, re- 
spectively. In order to monitor the Coulomb interaction 
effects, the scale-factor dependence has been studied by 
multiplying all the matrix elements Unn'mm' in eq.(73) 
with a factor A. Namely, the realistic value corresponds 
to A = 1. The effective model (71) with the above pa- 
rameters has been solved by PIRG. Technical details are 
found in Refs. 194,195. 

It has turned out that this compound shows very se- 
vere competitions as we see in Fig. 29. First, it lies on the 
verge of the Mott transition. Second, the ferromagnetic 
state is rather close in energy to the true ground state 
with the antiferromagnetic order. Third, spins and or- 
bitals order in a complicated pattern in the ground state 
as we see in Fig. 30 and candidates of the spin-orbital 
order are in severe competitions each other in the order 
of lOOK in energy. They have revealed rich orbital-spin 
physics arising from the competitions. 

Experimentally, transport and optical properties of 
this compound indicate either a very small Mott insulat- 
ing gap or semiconducting property with rapidly increas- 
ing resistivity with decreasing temperature^^^' in 
agreement with the above calculated results. The gap 
amplitude is nearly zero and it can easily be metallized 
by La doping. Recent experiments by dc susceptibil- 
ity and X-ray diffraction^^^ have suggested a transition 
around lOOK into a phase with antiferromagnetic and or- 
bital coupled order below this temperature, again in es- 
sential agreement with the above theoretical prediction. 
Recently, it has been proposed^^^ that the orbital-spin 
coupled order essentially described by the octupole or- 
der frequently discussed for /-electron systems^^^ might 
be stabilized when the spin-orbit coupling ignored in the 
available first-principles study are considered. 

In sharp contrast to this nearly insulating transport 
properties, the LDA calculation predicts a good metal- 
lic behavior (Fig. 28). On the other hand, the Hartree 
Fock approximation (HFA) predicts a clear ferromag- 
netic insulating phase at the realistic parameter values 
(see Fig. 29). LDA+U approach predicts results similar 
to HFA. The failure of single- Slater-determinant approx- 
imations as HFA and LDA+U is naturally understood 
because they relatively well describe a simple ferromag- 
netic state, while not the antiferromagnetic state with a 
nontrivial periodicity. Such a phase with large quantum 
fiuctuations can be described only by a more accurate 
solver such as PIRG beyond a single Slater determinant. 

All of the above agreement between the experiments 
and the present theory indicate that the approach of 
PIRG combined with the downfolding by using the LDA- 
GW scheme works well as a method for strongly corre- 
lated materials. From the viewpoint of the computational 
methods, Sr2V04 appears to offer a very severe and good 
benchmark for testing the accuracy in taking account of 
the correlation effects because of the severe competing 
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Fig. 28. (Color online) (Left panel) Electronic structure of 
Sr2V04 in LDA.194,195 (Right panel) Enlarged behavior of t2g 
bands computed from LMTO basis functions (solid light blue 
curves) and downfolded tight-binding bands (dot-dashed brown 
curves). The corresponding bands in the left panel are shown by 
arrows. The symbols denote the character of t2g bands in the 
F-point. The Fermi level is at zero energy. 




Fig. 29. Lowest energies per unit cell of total S=0 (open symbols) 
and ferromagnetic states (filled symbols) by quantum-number 
projected Hartree-Fock (triangles) and PIRG (circles) calcula- 
tions for the downfolded model of Sr2V04.-'^^^' ^'^^ (Inset): Lowest 
energies per unit cell of metallic (open squares) and insulating 
states (filled squares) by quantum- number projected PIRG for 
the downfolded model of Sr2V04. 

orders. 
5.2.5 YVOs 

YVO3 belongs to the family of transition-metal ox- 
ides with two valence electrons in the 3d orbitals {t2g 
manifold).^ The lattice structure is an orthorhombically 
distorted perovskite with the space group Phnm (four 
vanadium sites in a unit cell) at room temperatures. The 
GdFeOa-type distortion, rotation and tilting of the VOe 
octahedra are present, where the reduced V-O-V angle 
makes the narrow t2g bands. With lowering the temper- 
ature, it undergoes two successive phase transitions in 
both spin and orbital sectors. First, the G-type orbital 
ordering (00) appears at 200K with a structural change 
to the P2i/a symmetry, where a site with the dxy and dyz 
orbitals occupied and one with the the dxy and dzx are 
alternately arranged in three dimensions. The magnetic 
structure also shows the C-type spin ordering (SO) be- 
low 116K, where spins are aligned antiferromagnetically 




Fig. 30. (Color online) Ordered spin-orbital patterns in plane for 
Sr2V04 clarified by PIRG.-'^^'^' -"^^^ Ordered spin moment is pro- 
portional to the length of arrows. At each site, occupied or- 
bitals can be specified by a 3-dimensional unit vector in the ba- 
sis of t2g Wannier orbitals. Its xy^ yz, and zx components are 
given by (0.70,0.60,0.39), (0.51,0.80,0.31), (0.40,0.04,0.92), and 
(0.33,0.06,0.94), for the sites A, B, CI and C2, respectively 

in the a-b plane and ferromagnetically along the c-axis. 
With further lowering the temperature, the SO and 00 
simultaneously change at 77K, and the ground-state is 
the C-type 00 with the G-type SO.^^^'^^^ The crystal 
structure recovers the Phnm symmetry. In the charge 
sector, YVO3 is a typical Mott insulator with a large 
charge gap (~ leV). This is partly attributed to a large 
GdFeOs-type distortion, which reduces the bandwidth 
effectively. In addition, coupling to Jahn- Teller distor- 
tions is important in determining the orbital states. 

Electronic structure of YVO3 has been studied by the 
three-stage RMS scheme. The DFT-LDA calculations 
by using the local muffin-tin orbital basis has been ap- 
plied to derive the global band structure. The band struc- 
ture shown in Fig. 31 shows an isolated group of bands 
near the Fermi level mainly consisting of V 3d t2g atomic 
orbitals. The electron degrees of freedom far from the 
Fermi level are eliminated by a downfolding procedure 
leaving only the V 3d t2g Wannier bands as the low- 
energy degrees of freedom, for which a low-energy ef- 
fective model is constructed. This low-energy effective 
Hamiltonian is solved exactly by the PIRG method. 
It is shown that the ground state has the G-type spin 
and the C-type orbital ordering as we see in Fig. 32 in 
agreement with experimental indications. 

The indirect-charge gap is estimated to be 0.70 ± 0.07 
eV, which is smaller than the inferred experimental op- 
tical (direct) gap, but is consistent each other because 
the experiment has measured the direct gap while the 
gap in the calculation is the indirect gap. It has promi- 
nently improved the estimation compared to the LDA 
or GGA method. So far the indirect gap is not available 
experimentally. 

The LDA+PIRG results are thus all consistent with 
the available experimental results. In fact, this is the first 
result that reproduces an experimental charge gap as well 
as the spin and orbital ordering of YVO3 from the first- 
principles calculations. 

YVO3 and LaV03 have also been studied by the com- 
bination of the downfolding scheme and DMFT.^^^ It has 
been shown that the Jahn Teller and GdFe03 type dis- 
tortions are both crucial in reproducing the experimental 
orbital and magnetic orders at low temperatures. 
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Fig. 31. (Color online) Comparison of the band structure of 3d 
t2g orbitals computed from LMTO calculations (solid light blue 
lines) with the downfolded tight-binding model (dashed-dotted 
brown lines) for YVOs.^'-'^ 




Fig. 32. (Color online) Ordered spin- and orbital- patterns in the 
ground state of the PIRG solution for YVO3. The arrows rep- 
resent the magnetic local moment at each vanadium atom. The 
spins order antiferromagnetically in the G type, while the or- 
bitals order in the C type in agreement with the experiments. ^'-'^ 
The orbital states are shown in the form of the spatial electron 
distribution. 



5.2.6 Iron-based superconductors 

Recent discovery of iron based superconductors has 
renewed interest on high temperature superconductiv- 
ity/^' Several families of compounds are identified 
as superconducting materials, where Fe-3d conduction 
bands are commonly located near the Fermi level ac- 
cording to LDA band-structure calculations^^^'^^^ and 
their electrons are likely to form Cooper pairs. So far, 
the mechanism of superconductivity is not well un- 
derstood. In the family with ZrCuSiAs-type structure 
(called 1111 hereafter), SmFeAs(0,F) has the record of 
the highest superconducting critical temperature Tc ~ 56 
K (ref. 207) when fluorine is substituted with ~ 20% 
of oxygen as electron doping. There exist other fami- 
lies. BaFe2As2 with ThCr2Si2-type structure (called 122) 
shows the highest Tc ~ 38 K, when potassium is sub- 
stituted for ~ 40% of Ba as hole doping. LiFeAs 
and NaFeAs (called 111) with the PbFCl-type tetragonal 
structure indicate Tc - 18 K.^o^ 211 FeSe^^Tei-a, (called 
11) also shows superconductivity at Tc ~ lOK^^^'^^^) 
at ambient pressure and at Tc ~37 K under pressure (7 



GPa).2i4 

A common aspect of iron-based superconductors is the 
existence of antiferromagnetic order close to the super- 
conducting region except for the 111 family. However the 
ordered moment and pattern of the antiferromagnetism 
are strongly material dependent: LaFeAsO shows anti- 
ferromagnetic long-range order of the stripe type below 
Tat ~ 130 K with the Bragg point at (tt, 0) in the ex- 
tended Brillouin zone with a strongly reduced ordered 
moment ~0. 36-0. 63 /i^ as compared to the nominal sat- 
uration moment 4 /i^ for the high-spin 3d^ state. 2^^' 
Furthermore, LaFePO does not show an antiferromag- 
netic order and instead it undergoes a transition to the 
superconducting state at 4 K.^^'^ On the other hand, 
the 122-type (BaFe2As2) shows a relatively large ordered 
moment ~ 0.9 /i^ (refs. 218 and 219) and the 11-type 
(FeTe) indicates an even larger ordered moment ~ 2.0- 
2.25 /J.B at a different Bragg point, (7r/2, 7r/2).2i^'220 

Conventional LDA calculations of the 1111- 
type,22i 227 i22-type,228,229 iii_type,230 and 11-type 

compounds23^'232 g^ow a very similar band structure 
of the Fe 3d bands for all the compounds as we see in 
Fig. 33,^^^ where small electron pockets around M point 
and hole pockets around F point constitute semimetallic 
Fermi surfaces and the total widths of ten-fold Fe-3d 
bands are mostly around 4.5 eV. The local spin density 
approximation (LSDA) commonly predicted the antifer- 
romagnetic order for mother materials. 22^' 224, 226 rpj^^ 

stripe-type antiferromagnetic order is correctly repro- 
duced for the 1111-type. 224.226 However, the calculated 
ordered moment is unexpectedly too large (from 1.2 
to 2.6 ^g).223,224,226,233 -^^ contrast to much smaller 
ordered moment discussed above. The bicolinear order 
for FeTe is reproduced in the LSDA with the ordered 
moment ^ 2.25 /j,b) in agreement with the experimental 
results. 232 Diversity of the ordered moment ranging from 
zero to 2 /iB is surprising and not easily explained from 
the very similar band structure with semimetallic small 
pockets of the Fermi surface. Broad peak structures 
of magnetic Lindhard function calculated by using the 
LDA/GGA Fermi surface suggest severe competitions of 
different orderings.^^s, 225, 234 236 

Roles of electron correlations are not fully understood 
so far and are under debate. 2^'^ 242 Antiferromagnetic 
orders and fluctuations themselves revealed by the nu- 
clear magnetic resonance and other probes imply some 
electron correlation effects. 2^^' 2^3 S^iall fraction of the 
Drude weight244 247 -^^^ metallic behaviors^^' 2^^ 
support substantial electron correlation effects. Recent 
fluctuation exchange calculation suggests a substantial 
self-energy effect, where the validity of weak coupling 
and nesting picture becomes questionable. 2^^ 

Angle resolved photoemission spectroscopy2^^'25i has 
shown some correspondence to the LDA result of Singh 
et al.^^^ Fe-2p core-level spectra of X-ray photoemis- 
sion suggest rather itinerant character. 2^2, 253 However, 
some role of moderate electron correlations has also been 
claimed. ^^^'^^^ For FeSe, as we detail later, soft-Xray 
photoemission results2^^' 2^'^ appears to show a deviation 
from the LDA results and a crucial correlation effect. 2^^ 

In the superconducting phase, even the pairing sym- 
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parameters have been calculated by the ab initio down- 
folding and cRPA. So far, models for iron 3d orbitals 
{d model) and models including additional pnictogen or 
chalcogen p orbitals {dp or dpp models) have been de- 
rived. For the d model, the ratio of the averaged Hubbard 
diagonal onsite interaction U ^ 2.5 eV to a typical near- 
est neighbor transfer t ^ 0.3 eV in the downfolded model 
for LaFeAsO has been estimated to be U/t ^ 8-10 with 
the fivefold orbital degeneracy, indicating a moderately 
strong correlation. This moderately correlated nature has 
also been supported for the case of the 122-type, where 
U ^ 2.8 eV.^°^'^^^ For the case of the 11 compounds, the 
effective interaction is even larger as C7 ~ 4.2 eV for FeSe 
and 3.4 eV for FeTe. In Fig. 35, comparisons of the derived 
ab initio model parameters are shown. We note that the 
effective Coulomb interaction for the low-energy down- 
folded model estimated in these works is different from 
the interaction observed by experimental probes such as 
the X-ray photoemission.^^^'^^^ The measured interac- 
tion parameters are resulted from the further screening 
by the 3d electrons excluded in the model construction. 



Fig. 33. Electronic band structures of six iron-based supercon- 
ductors obtained by DFT-LDA.^o^ The Ki - K5 points in 
BaFe2As2 are Ki = ^(i,0,0),X2 = ^(ii0).^3 = 
^(0,0,_^),K4= ^(iO,^),K5 = ^(ii^), respectively. 
Energy is measured from the Fermi level. 
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FeSe 



FeTe 





metry itself is highly controversial and no consensus 
has been reached. Although nodeless superconductivity 
is suggested, ^^^^^^ temperature dependence of nuclear- 
magnetic-relaxation time Ti below Tc roughly scaled by 
without the Hebel-Slichter peak implies unconven- 
tional superconductivity driven by nontrivial electron- 
correlation effects. For example, orbital dependent 
gaps with sign-changing and fully-gapped s± symmetry 
has been proposed. The gradual suppression of the su- 
perconducting transition temperature by Co doping into 
the Fe site was reported to be explained by the 5-wave 
singlet pairing without the sign change. ^^^^^^ Although 
overall experimental results suggest noticeable correla- 
tion effects, realistic roles on the pairing are not well 
established and controversial. 

To understand properties and mechanisms of mag- 
netism and superconductivity in iron based supercon- 
ductors, and to distinguish what are common and what 
are family dependent, effective low-energy models of 
these families have been derived^^^' from first prin- 
ciples along the line of RMS based on the three-stage 
scheme. In the procedure, the LDA band struc- 
ture was calculated as we see in Fig. 33 and the maximally 
localized Wannier functions are constructed as we see ex- 
amples in Fig. 34, from which the transfer and interaction 



Fig. 34. (Color online) Isosurface of maximally localized Wan- 
nier function at ±0.02 a.u. for Fe — orbital in d model of 
LaFeAsO (left), FeSe (middle), and FeTe (right ).206 This illus- 
trates how the Wannier spread shrinks from LaFeAsO to FeTe. 
The dark shaded surfaces (color in blue) indicate the positive 
isosurface at +0.02 and the light shaded surfaces (color in red) 
indicate —0.02. 

The systematic change in the model parameters is un- 
derstood from the structural differences. A key quantity 
for understanding the systematic evolution from the 1111 
to the 11 families is the height of the pnictogen/chalcogen 
layer h measured from the iron plane. The height h 
increases from LaFePO (1.13 A), LaFeAsO (1.32 A), 
BaFe2As2 (1.36 A) FeSe (1.47A) to FeTe (1.77A). It was 
pointed out in the early stage that the electronic band 
structure is altered significantly by changing Z^. 222, 266 
There is also a work claiming that the spin and charge 
susceptibilities are sensitive to /i.^^^ In terms of the corre- 
lation strength, h controls the spatial extent of the Wan- 
nier orbitals and the strength of screening effect. The 
smaller h enhances the hybridization of the Fe 3d or- 
bitals with the pnictogen/chalcogen p orbitals leading to 
extended Wannier orbitals as is shown in Fig. 34. This 
makes the bare onsite interaction small. The smaller h 
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Fig. 35. (Color online) Material dependence of parameters for models for iron 3d orbitals.206 The average of the onsite intraorbital 
effective Coulomb interactions (f7), the average of the offsite effective Coulomb interactions between the neighboring Fe sites at the 
same orbitals(1/), the average of the onsite effective exchange interactions (J), the maximum value of the transfer integrals between the 
neighboring Fe sites [tn = tii(l/2, —1/2, 0)] and between the next-nearest neighbor [t^^ = ^44(1, 0,0)], U /t, and ^44/^11 are compared. 
The subscripts of tn and ^44 are orbital indices; 1 for xy and 4 for zx. t is the orbital average of the largest nearest d-d transfer 
integrals. 



also makes the pnictogen/chalcogen p level closer to the 
Fermi level, which enhances the screening of the Coulomb 
interaction by the p bands. Furthermore, the number of 
p bands contributing to the screening decreases in sim- 
pler compounds such as the 11 compounds, which re- 
duces the screening channels and enhance the effective 
interaction for large h. The difference is similar in the 
effective models containing p orbitals of As, Se or Te {dp 
or dpp model), where U ranges from ~ 4 eV for the 1111 
family to ~ 7 eV for the 11 family. The exchange inter- 
action J has a similar tendency. The family dependence 
of models indicates a wide variation ranging from weak 
correlation regime (LaFePO) to substantially strong cor- 
relation regime (FeSe). This variety of the electron cor- 
relation brings about the diversity of physical properties 
observed in different families of iron based superconduc- 
tors in spite of similar band structures. 

In fact it has been pointed out that FeSe may show 
the Hubbard splitting of the bands as a clear indication 
of the strong correlation effects. The available exper- 
imental results are consistently analyzed from this per- 
spective. For the 1111 family, the correlation effects have 

been analyzed in more detail as a moderately correlated 
system.239,267,268 

The larger h also explains why the ten-fold 3d bands 
for the 11 family are more entangled with the smearing 



of the "pseudogap" structure above the Fermi level ob- 
served in the 1111 family. While the family-dependent 
semimetallic splitting of the bands primarily consists 
of dyz/dzx and dx2_y2 orbitals, the size of the pseudo- 
gap structure is controlled by the hybridization between 
these orbitals and d^y j d-^^'^ A large hybridization in 
the 1111 family generates a large "band- insulating" -like 
pseudogap {hybridization gap), whereas a large h in the 
11 family weakens them, resulting in a "half- filled" like 
bands of orbitals. This enhances strong correlation effects 
in analogy with Mott physics and causes the orbital se- 
lective crossover in the three orbitals. On the other hand, 
the geometrical frustration f /t, inferred from the ratio of 
the next-nearest transfer t' to the nearest one t of the d 
model is relatively larger for the 1111 family than FeTe. 

By using the ab initio model, magnetic properties have 
been analyzed by low-energy solver s.^^ The magnetic 
transition with the correct pattern has been reproduced 
by the many- variable variational Monte Carlo calcula- 
tions with a quantitative agreement of the ordered mo- 
ment. In case of LaFeAsO, the unusually small moment 
has been naturally understood by the proximity to a 
quantum critical point between a paramagnetic metal 
and an antiferromagnetic metal. VMC results are shown 
in Fig. 36 for the ordered magnetic moment as a func- 
tion of the scaled interaction strength to monitor the 
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interaction effects. The ordered moment shows system- 
atic evolution from the quantum critical point near the 
interaction corresponding to LaFaAsO and is in quanti- 
tative agreement with the experimental results shown by 
crosses without an adjustable parameter. Here, the ab 
initio model was constructed for LaFeAsO correspond- 
ing to A = 1 and the parameter A for other compounds 
is determined by the ratio of the averaged onsite intraor- 
bital interaction between the compound and the refer- 
ence system LaFeAsO. This ratio was calculated from 
the ab initio model parameters obtained in Ref.206. The 
parameter A is obtained from the original scaling param- 
eter A by considering the La 4/ screening ignored in the 
model construction by Nakamura et al.^^^ and by consid- 
ering the interlayer screening effects discussed in §3.8.^^ 
It is remarkable that the ab initio downfolded models for 
various different families of the iron superconductors are 
all within a few percent of errors given by a single effec- 
tive Hamiltonian with a single parameter A. The robust 
metallic behavior in a large interaction region {U /t ~ 10) 
is also understood from the existence of two Dirac cones 
in the dispersion near the Fermi level. The metal is pro- 
tected as long as the Dirac cones are retained. The Dirac 
cones can be annihilated in pair only by a large magnetic 
moment (> S/ie)- 
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Fig. 36. (color online): Magnetic ordered moment m(gpeak) as a 
function of the interaction parameter A scaling the overall in- 
teraction amplitude calculated by VMC and extrapolated to the 
thermodynamic limit for the ab intitio model of LaFeAsO (open 
circles). @ The effective scaling factor A is obtained by fur- 
ther considering the interlayer screening effect (see the text). 
Experimentally observed materials dependence at correspond- 
ing A is also shown by crosses for LaFePO, ^-'^'^ LaFeAsO, ^-"^^ 
BaFe2Se2^^^'2^^ and FeTe.. ^i^, 220 Quantum critical point of the 
AF transition appears at slightly below A = 0.75(A ~ 1). 



5.3 Organic conductor 

Families of organic conductors provide us with another 
type of strongly correlated electron systems. Usually, 
a unit cell of molecular crystals contains many atoms. 
However, in many cases, a unit cell contains only a small 
number of molecules and structures are simple in terms 
of molecular stackings. Electrons in a single molecule 
occupy molecular orbit als, which has normally spread 
over the molecule. Such molecular orbitals have a small 
overlap with those on neighboring molecules, when the 
molecules are stacked to form the bulk crystal. Typi- 
cal bands have simple structures near the Fermi level, 
where they consist only of lowest unoccupied molecu- 
lar orbital (LUMO)) and highest occupied molecular or- 
bital (HOMO). The HOMO and LUMO bands are in 
many cases isolated from other bands. This makes the 
effective Coulomb interaction between electrons on the 
HOMO and LUMO orbitals poorly screened by other 
bands. The small overlap of the molecular orbitals be- 
tween the neighboring molecules makes the HOMO and 
LUMO bandwidths small. These two factors contribute 
to make the organic conductors mostly strongly corre- 
lated electron systems. 

5.3.1 n-ET conductor 

ET-type conductors are synthesized as a family of such 
organic conductors. (ET)2X with a number of choices 
of anions X, offer a variety of prototypical behaviors of 
correlated electron systems with two-dimensional (2D) 
anisotropics. Examples range from correlated metals 
with superconductivity at low temperatures and com- 
peting charge orderings to Mott insulators either with a 
quantum spin liquid or with antiferromagnetic, charge- 
density or spin-Peierls orders. Intriguing Mott transi- 
tions are also found. They have all been studied exten- 
sively at a front of research for unconventional quantum 
phases and quantum critical phenomena. Among them, 
/^-type ET conductors have stacking of dimerized pair of 
ET molecules. The dimerization causes splittings of the 
HOMO and LUMO bands into bonding and antibond- 
ing bands. Since holes are quarter filled (electrons are 
three-quarter filled) at the HOMO band for (ET)2X, af- 
ter the dimerization splitting, the Fermi level is normally 
located at the antibonding HOMO band at half filled. 

In particular, an unconventional nonmagnetic Mott- 
insulating phase is found near the Mott transition in the 
n-iype structure of ET molecules, /t:-(ET)2Cu2(CN)3 re- 
ferred to as /t:-CN. Although this compound is a Mott in- 
sulator, no magnetic order is identified down to the tem- 
perature T=0.03 K, four orders of magnitude lower than 
the antiferromagnetic spin-exchange interaction J~250 

270 rpj^^ emergence of the quantum spin liquid near 
the Mott transition has been predicted in earlier numer- 
ical studies,^^'^^'^^^ while the full understanding of the 
spin liquid needs more thorough studies. It is also cru- 
cially important to elucidate the real relevance of the 
theoretical findings to the real /^-ET compounds. Most 
of numerical^^^' and theoretical^^^ studies have also 
been performed for a simplified single-band 2D Hubbard 
model based on an empirical estimate of parameters com- 
bined with extended Hiickel calculations. ^^^'^^^ A more 
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realistic description of k,-ET compounds was certainly 
needed beyond the empirical model. 

Another fundamental finding achieved in this series of 
compound is the unconventional Mott transition found 
for X=Cu[N(CN)2]Cl under pressure. The novel uni- 
versality class of the Mott transition is in good agreement 
with the theoretically revealed marginal quantum criti- 
cality at the meeting point of the symmetry breaking and 
topological change. ^^''^^^ Because of its significance to 
the basic understanding on the physics of quantum criti- 
cality, the relevance of theoretical concept to the experi- 
mental observation needs to be further examined on real- 
istic and first-principles grounds. Furthermore, an uncon- 
ventional superconductivity is found in the metallic sides 
of these compounds at low temperatures (T<Tc~ 10- 
13K), where the mechanism is not clearly understood 
281,282 xhese outstanding properties of a^-ET com- 
pounds have urged systematic studies based on realistic 
basis. As mentioned above, however, the first-principle 
studies are Hmited^^^ and most of the studies so far were 
performed using the empirical models inferred from the 
Hiickel studies. 

The three-stage RMS method has recently been 
applied to tz-ET conductors. Effective low-energy 
models have been derived for two contrasting com- 
pounds, spin-liquid /^-CN and superconducting com- 
pound X=Cu(NCS)2 abbreviated as /^-NCS,^^^ to get 
insights into the whole series of f<i-{ET)2X compounds 
from metals to Mott insulators. 

The global band structures obtained by GGA are 
shown in Figs. 37(a) and (b). They clearly and commonly 
show that the antibonding HOMO bands are isolated 
near the Fermi level as is anticipated. Then the down- 
folding to an effective single-band model for the HOMO 
antibonding band has been performed after construct- 
ing the maximally localized Wannier orbital shown in 
Fig.37(c).284 

The parameters of the downfolded model are listed 
in Table H for a^-NCS and a^-CN with the notation 
of the transfer in Fig. 37(d). It contains dispersions of 
the highest occupied Wannier-type molecular orbitals 
with the nearest neighbor transfer t~0.067 eV for a 
metal X=Cu(NCS)2 and 0.055 eV for a Mott insulator 
X=Cu2(CN)3, as well as the onsite screened Coulomb 
interactions. It shows a substantial difference from the 
previous simple extended Hiickel results:^^^'^^^'^^^ The 
derived parameters indicate that (i) the geometrical frus- 
tration parameter \t' /t\ is substantially smaller than the 
extended Hiickel results and a^-CN estimated at \t' /t\^O.S 
has turned out to be away from the right triangular struc- 
ture^^^ and (ii) the onsite Coulomb repulsion {U^O.S eV 
characterized by /7/t~12-15) is unexpectedly large com- 
pared to the Hiickel estimate given by U/t^7-S^ while 
the intersite Coulomb interaction was found to be also 
appreciable as we see in Figs. 38(a) and (b). 

6. Concluding Remark and Outlook 

We have reviewed recent rapid advance in understand- 
ing electronic structures of strongly correlated electron 
systems by utilizing the electronic hierarchical structure. 
The method starts from obtaining the global electronic 
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Fig. 37. (Color) Ab initio GGA band structures (red line) of 
As:-(BEDT-TTF)2Cu(NCS)2 (a) and At-(BEDT-TTF)2Cu2(CN)3 
(b).^^^ The zero of energy corresponds to the Fermi level. The 
blue dotted dispersions are obtained by the four transfer pa- 
rameters listed in TABLE I. (c) Maximally localized Wannier 
function of as:-(BEDT-TTF)2Cu(NCS)2 constructed for effective 
low-energy model. The amplitudes of the contour surface are 
+1.5/v^ (blue) and — 1.5/v^ (red), where v is the volume of the 
primitive cell. S, C, H, N, and Cu nuclei are illustrated by green, 
yellow, silver, blue, and red spheres, respectively, (d) Schematic 
network of transfers in the triangular lattice. 



Table II. List of the parameters in a form of the single-band ex- 
tended Hubbard Hamiltonian in eq. (71) for ^c-(ET)2X. 

X=Cu(NCS)2 X=Cu2(CN)3 



ta (meV) -64.8 -54.5 
h (meV) -69.3 -54.7 
tc (meV) 44.1 44.1 
td (meV) -11.5 - 6.8 
U (eV) 083 0.85 



structure by DFT or GW procedure. Then low-energy 
effective models are derived by the downfolding elimi- 
nating the degrees of freedom away from the Fermi level. 
This is achieved first by extracting localized Wannier or- 
bitals and constructing ab initio tight binding models in 
real space. The interaction parameters are obtained by 
counting the screening with the constrained RPA. The 
resultant models are solved by highly accurate solvers 
such as the various Monte Carlo methods, path-integral 
renormalization group and dynamical mean field theory. 
This scheme opens a way to understand electron correla- 
tions even when the single-particle picture breaks down. 

The applications to transition metal compounds in- 
cluding iron-based superconductors and perovskite ox- 
ides as well as to organic conductors have already proven 
its quantitative accuracy without any adjustable parame- 
ters and shown that a new powerful method has emerged. 
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Fig. 38. (Color online) Screened Coulomb interactions of k- 
(BEDT-TTF)2Cu(NCS)2 (a) and as:-(BEDT-TTF)2Cu2(CN)3 
(b) as a function of distance between two centers of maximally 
localized Wannier orbitals, calculated by cRPA illustrated by 
circles (red). Crosses (green) show the bare interactions. -^^^ The 
solid and dotted curves denote 1/r and l/(Ar), with a fitting 
parameter A ~ 5 hartree^^bohr"-*^ in both of the compounds. 



Severely corapeting orders as well as quantura and raany- 
body fluctuations are now under controlled treatment by 
this approach beyond the mean-field and one-body pic- 
tures, while it is still on the way toward further growth 
with diverse potentiality of improvement. Recent ap- 
plications to alkali-cluster- loaded zeolites^^^' show a 
wide feasibility and potential of the methods including 
the compounds with a large unit cell. Another direction 
of the application is surfaces and interfaces, which do not 
retain the bulk translational symmetry. In such complex 
systems, the number of bands near the Fermi level can 
be very large because of large unit cells. 

When electron correlation effects are so large that the 
degrees of freedom far from the Fermi level are modified, 
one has either to include such a part under correlation ef- 
fects into the low-energy effective model, or has to solve 
them selfconsistently with a feedback to a high-energy 
downfolded part as is illustrated as the broken arrow in 
Fig.l. Although preliminary results have been reported 
for LDA+DMFT and GW+DMFT, applications are so 
far limited because of computationally demanding iter- 
ations. The total selfconsistency is certainly a future di- 
rection of challenge. 

We have considered only the electronic degrees of free- 
dom and the atomic structure is assumed to be given 
in this article. A future important subject is to com- 
bine with the structural optimization for the goal of the 
real first principles scheme. In addition, phonon de- 
grees of freedom are in general coupled in a low-energy 
scale^^^^^^ and it determines a number of important 
properties including phonon mediated superconductiv- 
ity, ferroelectricity and charge ordering. 

Another intriguing problem is dynamical processes 
far from equilibrium and relaxation phenomena. Pho- 
toinduced transitions are typical examples of future is- 
sue. Experimentally, time resolved photoemission will 
open powerful probes and tools for new situation of 
nonwquilibrium phenomena. Depending on the energy 
and time scales, we need to develop more involved but 
tractable framework with low-energy solvers for nonequi- 
librium. 
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